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)})