This paper is concerned with a high-order numerical scheme for nonlinear systems of second-order boundary value problems (BVPs). First, by utilizing quasi-Newton’s method (QNM), the nonlinear system can be transformed into linear ones. Based on the standard Lobatto orthogonal polynomials, we introduce a high-order Lobatto reproducing kernel method (LRKM) to solve these linear equations. Numerical experiments are performed to investigate the reliability and efficiency of the presented method.