Cohesive zone models are being increasingly used to simulate fracture and debonding processes in metallic, polymeric, and ceramic materials and their composites. The crack initiation process as well as its actual stress and damage distribution beyond crack tip are important for understanding fracture of materials and debonding of adhesively bonded joints. In the current model, a natural boundary condition based method is proposed, and thus the concept of extended crack length (characteristic length ) is no longer required and more realistic and natural local deformation beyond crack tip can be obtained. The new analytical approach, which can consider both crack initiation and propagation as well as local deformation and interfacial stress distribution, can be explicitly obtained as a function of the remote peel load with the given bilinear cohesive laws. An intrinsic geometric constraint condition is then used to solve the remote peel load . The nonlinear response in both the ascending and descending stages of loading is accurately predicted by the current model, as evidenced by a comparison with both experimental results and finite element analysis results. It is found that the local deformation and interfacial stress beyond crack tip are relatively stable during crack propagation. It is also found that, when the cohesive strength is low, it has a significant effect on the critical peel load and loadline deflection. In principle, the approach developed in the current study can be extended to multilinear cohesive laws, although only bilinear cohesive law is presented in this work as an example.