def rk4_step(y,t,h,deriv):
    """
    General RK4 stepper for
    N coupled differential eq's

    Input:
     - y:      Array containing the y(t)
     - t:      Which time are we talking about?
     - h:      Stepsize
     - deriv:  Function that returns an array
               containing dy/dt for each y at time t,
               and takes as arguments an y-array, and time t.
    """
    
    k1 = h*deriv(y,t);
    k2 = h*deriv(y+k1/2.0,t+h/2.0)
    k3 = h*deriv(y+k2/2.0,t+h/2.0)
    k4 = h*deriv(y+k3,t+h)

    return y + (k1 + 2*(k2+k3) + k4)/6.0