In this study, a finite element strategy is implemented for more accurate simulation of nonlinear failure behavior in multidirectional laminated composites. The considered approach is based on the parameters which could affect the SERR2. Accordingly, interaction between delamination and matrix cracking as well as mixed stress distribution around the fiber due to mismatch in properties of adjacent plies with different fiber angle orientation influence the strain energy release rate. Thus, in the present research the prediction of the damage behavior of cross-ply laminates is improved considering the effects of matrix cracking and inhomogeneous fiber-matrix modeling. For simulation of the delamination and the matrix cracking damage, cohesive interface elements are implemented. They have been used perpendicularly inside the 90o plies for the matrix cracking. CZM3 approach, ABAQUS subroutine and Python programming were utilized for automatic modeling of matrix cracking. The resultant load–displacement curve in DCB4 test of cross-ply laminates indicated that by considering delamination and matrix cracking simultaneously, load–displacement curve of cross-ply laminates became nonlinear earlier and the obtained results had an acceptable compatibility with experimental results. The obtained results for double cantilever beam (DCB) test specimen with cross-ply lay-up indicate that the numerically reproduced GIc parameter of 90/0 interface can be improved about 20% using this simulation. Moreover, by simulating periodic matrix cracking inside 90o plies, crack jumping, stick-slip phenomenon and local hardening behavior are anticipated such as experiments.