The debonding process in FRCM (Fiber Reinforced Cementitious Matrix) strengthened brittle substrates, such as concrete and masonry, is characterized by complex failure modes that include debonding at both the external and internal interfaces of the composite, as well as damage to the mortar matrix. This study addresses this problem classically as a mode-II fracture process by accounting for the one-dimensional non-linear interfacial shear stress-slip relationship combined with matrix axial non-linearity. The approach employed considers the longitudinal displacements of both the inner and outer layers of mortar, along with that of the fiber textile, as unknowns, and analytically solves the second-order ordinary differential equation (ODE) system that governs the debonding problem. Then, the boundary value problem (BVP) is transformed into an initial value problem (IVP) and a double shooting method is proposed to find displacements and stresses along the bonded length in the presence of non-linearity. The appropriate initial values of matrix layer displacements, satisfying the required boundary conditions, are efficiently determined through a two-dimensional bisection method working on rectangular and triangular patches. The softening behavior of both mortar layers and interfaces between mortar and fiber are taken into account using jagged constitutive relationships to enhance convergence. The numerical calculation speed, robustness, and key simulation parameters are thoroughly discussed. Moreover, the proposed numerical approach is compared with existing experimental data and models, on both concrete and masonry FRCM-reinforced specimens. The influence of mortar failure and interfacial strength is also explored. The model demonstrates its ability to effectively reproduce the global bond behavior observed in experimental studies, while capturing the local behavior and simulating various failure modes observed in the tests.