Methods to solve systems of first-order differential equations are generalizations of the methods for a single first-order equation. For example, using the notation that \(\omega_{ij}\) denotes the approximation at \(t=t_j\) of the \(i^{th}\) solution \(y_i(t)\) for each \(j=0 \dots N\) and \(i= \dots m\), then the fourth order Runge-Kutta method may be written as:

\begin{equation} \begin{aligned} \omega_{1,0} & = \alpha_1 \quad \omega_{2,0} = \alpha_2 \dots \omega_{m,0} = \alpha_m \\ k_{1,i} & = hf_i(t_j, \omega_{1,j}, \dots, \omega_{m,j}) \\ k_{2,i} & = hf_i \left(t_j + \frac{h}{2}, \omega_{1,j} + \frac{1}{2}k_{1,1}, \omega_{2,j} + \frac{1}{2}k_{1,2}, \dots, \omega_{m,j} + \frac{1}{2}k_{1,m} \right) \\ k_{3,i} & = hf_i \left(t_j + \frac{h}{2}, \omega_{1,j} + \frac{1}{2}k_{2,1}, \omega_{2,j} + \frac{1}{2}k_{2,2}, \dots, \omega_{m,j} + \frac{1}{2}k_{2,m} \right) \\ k_{4,i} & = hf_i \left(t_j + h, \omega_{1,j} + k_{3,1}, \omega_{2,j} + k_{3,2}, \dots, \omega_{m,j} + k_{3,m} \right) \\ \omega_{i,j+1} & = \omega_{i,j} + \frac{1}{6}(k_{1,i} + 2k_{2,1} + 2k_{3,i} + k_{4,i}) \\ \end{aligned} \label{eq:sys1} \end{equation}for each \(i=1 \dots m\). Note that all the values \(k_{1,1}, k_{1,2}, \dots, k_{1,m}\) must be computed before any of the terms of the form \(k_{2,i}\) can be determined. In general, each \(k_{l,1}, k_{l,2}, \dots, k_{l,m}\) must be computed before any of the expressions \(k_{l+1,i}\). The other one-step methods can be extended to systems in a similar way. When error control methods are extended, each component of the numerical solution \(\omega_{1,j}, \omega_{2,j}, \dots, \omega_{m,j}\) must be examined for accuracy. If any of the components fail to be sufficiently accurate, the entire numerical solution must be recomputed. Multistep methods may also be extended to systems of equations.