We provide of a method to integrate first order non-linear systems of
differential equations with variable coefficients. It determines
approximate solutions given initial or boundary conditions or even for
Sturm-Liouville problems. This method is a mixture between an iterative
process, a la Picard, plus a segmentary integration, which gives
explicit approximate solutions in terms of trigonometric functions and
polynomials. The segmentary part is particularly important if the
integration interval is large. This procedure provide a new tool so as
to obtain approximate solutions of systems of interest in the analysis
of chemical reactions. We test the method on some classical equations
like Mathieu, Duffing quintic equation or Bratu's equation and have
applied it on some models of chemical reactions.