In this paper, an optimization algorithm is presented for solving systems of s ingular 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 u sing continuous genetic algorithm in the sense that each of the derivatives is replaced by an appropriate difference quotient appro ximation. Two numerical experiments are carried out to verify the mathematical results, and the theoretical statements for the solutions are su pported by the results of numerical experiments. Meanwhile, the statistical analysis is provided in order to capture the behavior of the solu tions 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.