A fast bisection procedure embedded in a standard ODE45 solver is proposed for three different kinds of FRCM experimental tests, namely (1) coupon tensile test, (2) pull-out test and (3) single-lap shear test. In order to predict global and local behavior of FRCM in such kind of experiments, a mono-axial differential equations model can be derived, requiring in the most general case the numerical solution of Boundary Value Problems BVPs, involving as independent variables mortar and fiber displacement fields and their first derivatives. Assuming a nonlinear behavior of both matrix and reinforcement-matrix interface, implicit and general purpose BVPs exhibit poor numerical stability and slow convergence. Furthermore, the initial guess of the solution is an integral part of solving a BVP, making the approach even less appealing, especially to reproduce FRCM behavior in the whole possible range of deformations, where severe softening and snap back phenomena are encountered. In the paper, the BVP is transformed in an Ordinary Differential Equation ODE system with initial values, i.e. into an Initial Value Problem IVP, by means of a classic shooting method where guess values of the matrix displacement on the opposite side of the loaded edge are assumed and varied until the desired boundary condition on the loaded edge are obtained. A classic bisection procedure allows to converge quickly to the exact value of matrix displacement to assign in order to recover quickly the boundary condition of the original BVP problem imposed on the opposite side. The advantage -common to all shooting methods available- is that a BVP is transformed in a sequence of ODE problems with initial conditions, which can be solved by means of explicit Runge Kutta kernels. Any kind of non-linearity can be assumed for matrix –for which a particular saw tooth approximation of the real behavior in tension is proposed to strengthen further convergence- and matrix-reinforcement interface, for which both a newly proposed exponential law and a stepped softening relationship are adopted. The numerical procedure is tested on the aforementioned three different experimental setups, showing excellent fitting capability of both the experimental behavior and previously presented theoretical models, where available.