Updated 2012-05-08: There was a pretty serious mistake in my original implementation: it didn't work for more than one dependent variable.

This is a Python implementation of RK4. It can be used to solve systems of first order ordinary differential equations (and, by extension, higher order ordinary differential equations).

For an implementation specialized for second order ordinary differential equations, see my older post.

Here is the main function:

    def rk4(x, h, y, f):
        k1 = h * f(x, y)
        k2 = h * f(x + 0.5*h, y + 0.5*k1)
        k3 = h * f(x + 0.5*h, y + 0.5*k2)
        k4 = h * f(x + h, y + k3)
        return x + h, y + (k1 + 2*(k2 + k3) + k4)/6.0
    

"That's all?" you may ask. Yes, it is.

Here is an example, using numpy for numerical arrays. Note that numpy overloads the + and * operations on arrays so that we don't have to iterate over each element by hand.

    import numpy

    def damped_spring(t, state):
        pos, vel = state
        stiffness = 1
        damping = 0.05
        return numpy.array([vel, -stiffness*pos - damping*vel])

    t = 0
    dt = 1.0/40
    state = numpy.array([5, 0])
    print('%10f %10f' % (t, state[0]))

    while t < 100:
        t, state = rk4(t, dt, state, damped_spring)
        print('%10f %10f' % (t, state[0]))
    

If you run this program, it will output a line for each timestep indicating the time that has passed and the current position of the spring. If you plot it, you'll get something like this:

    sh$ python example.py > output.dat
    sh$ gnuplot
    gnuplot> plot 'output.dat' using 1:2 with lines