The central piece of Rungu-Kutta method is the approximation of the increasement of the function. In 1st order ODE,

$\dot{y} = f(x,y), \dot{y}(x_0) = y_0$

In a special case of Rungu-Kutta of order 4 (RK4), there are 2 array $b_i$ and $c_i$, so that

$dx = h, dy = h \sum\limits_{i=1}^{4} b_i dy_i$

$dy_i = f(x+ c_i dx, y+c_i dy_{i-1})$

where $c_i$ is ranging from 0 to 1, $\sum\limits_{i=1}^{4} b_i =1$.

In RK4,

$c = (0, \frac{1}{2}, \frac{1}{2}, 1)$

$b = \frac{1}{6}(1,2,2,1)$

The $c_1 = 0$ is a must, otherwise, we have to define $dy_0$. There should be some methods to obtain an optimum values for $c$ and $b$, but I don’t know.

For 2nd order ODE

$\frac{d^2y}{dx^2} + a \frac{dy}{dx} = f(x,y), y(x_0) = y_0, \frac{dy}{dx}|_{x_0} = z_0$

change to

$\frac{dy}{dx} = z, \frac{dz}{dx}=F(x,y,z)$

These equation are the similar 1st order ODEs.

$dx = h, dy = \sum\limits_{i=1}^{4} b_i dy_i, dz = \sum\limits_{i=1}^{4} b_i dz_i$

$dy_i = h( z+ c_i dz_{i-1})$, where is the z for $y_n$ step.

$dz_i = h F(x+ c_i dx, y+c_i dy_{i-1}, z+c_i dz_{i-1})$

The $dy_i$ is using $dz_i$

$x_{n+1} = x_n + dx$

$y_{n+1} = y_n + dy$

$z_{n+1} = z_n + dz$

This can be generalized to any order ODE by decoupling the ODE into

$\frac{dy}{dx} = x,..., \frac{du}{dx}=F(x,y,...,u)$

the equation for $du_i$ is

$du_i = h F(x+ c_i dx, y+c_i dy_{i-1}, z+c_i dz_{i-1}, ...., u+ c_i du_{i-1})$

And for all the intermediate variable $y, z, w,.... = O^{(k)}$

$dO_i^{(k+1)} = h( O^{(k)} + c_i dO_{i-1}^{(k)})$