In this paper, an optimization algorithm is presented for solving systems of singular boundary value problems. In this technique, the system is formulated as an optimization problem by the direct minimization of the overall individual residual error subject to the given constraints boundary conditions, and is then solved using continuous genetic algorithm in the sense that each of the derivatives is replaced by an appropriate difference quotient approximation. Two numerical experiments are carried out to verify the mathematical results, and the theoretical statements for the solutions are supported by the results of numerical experiments.Meanwhile, the statistical analysis is provided in order to capture the behavior of the solutions and to discover the effect of system parameters on the convergence speed of the algorithm. The numerical results demonstrate that the algorithm is quite accurate and efficient for solving systems of singular boundary value problems. © 2014 NSP Natural Sciences Publishing Cor.