I was trying to iterate a system of differential equations and decided to implement it directly in Python to refresh my memory. This is one way to find the approximation of solutions of ordinary differential equations, give some initial value.

Higher-order methods (such as RK4) are generally used due to their superior convergence and stability properties as well as larger step size.

Given a differential equation with an initial condition  $\dot{y} (t) = f(t,y)$ with $ y(t_0) = y_0$ where we would like to find $y$. We pick a stepsize $h$ and update the value of $y(t+h)=y_{n+1}$ using

$$ y_{n+1} = y_n + \frac{h}{6} (k_1 + 2k_2 + 2 k_3 + k_4 )$$

$$t_{n+1}=t_n + h$$

where the k functions are

$$k_1 = f(t_n, y_n) $$

$$k_2 = f(t_n+\frac{1}{2} h, y_n+\frac{h}{2} k_1) $$

$$k_3 = f(t_n+\frac{1}{2} h, y_n+\frac{h}{2} k_2) $$

$$k_4 = f(t_n+ h, y_n+h k_3) $$

This can be extended to solve several functions simultaneously of the kind

$$ \dot{y} = f(t,y,z) \quad , \quad \quad \dot{z} = g(t,y,z) $$

The RK4 then becomes

$$ y_{n+1} = y_n + \frac{h}{6} (k_1 + 2k_2 + 2 k_3 + k_4 )$$

$$ z_{n+1} = z_n + \frac{h}{6} (l_1 + 2l_2 + 2 l_3 + l_4 )$$

$$t_{n+1}=t_n + h$$

where the k functions are

$$k_1 = f(t_n, y_n, z_n) $$

$$l_1 = g(t_n, y_n, z_n) $$

$$k_2 = f(t_n+\frac{1}{2} h, y_n+\frac{h}{2} k_1, z_n+\frac{h}{2} l_1) $$

$$l_2 = g(t_n+\frac{1}{2} h, y_n+\frac{h}{2} k_1, z_n+\frac{h}{2} l_1) $$

$$k_3 = f(t_n+\frac{1}{2} h, y_n+\frac{h}{2} k_2, z_n +\frac{h}{2} l_2 ) $$

$$l_3 = g(t_n+\frac{1}{2} h, y_n+\frac{h}{2} k_2, z_n +\frac{h}{2} l_2 ) $$

$$k_4 = f(t_n+ h, y_n+h k_3, z_n+h l_3) $$

$$l_4 = g(t_n+ h, y_n+h k_3, z_n+h l_3) $$

Extension to more than one simultaneous differential equation from JBM.

Neat implementation using lambda functions on Rosseta Code.

def RK4(f):
    return lambda t, y, dt: (
        lambda dy1: (
        lambda dy2: (
        lambda dy3: (
        lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
        )( dt * f( t + dt , y + dy3 ) )
        )( dt * f( t + dt/2, y + dy2/2 ) )
        )( dt * f( t + dt/2, y + dy1/2 ) )
        )( dt * f( t , y ) )

You can find the whole code here.

If you are interested using this for something useful, then take a look at the RKF45 Runge-Kutta-Fehlberg ODE Solver written in Python by Peter Monk. It is also available in Fortran, C and Matlab.