def rk4(derivs, y0, t):
try: Ny = len(y0)<EOL>except TypeError:<EOL><INDENT>yout = np.zeros( (len(t),), np.float_)<EOL><DEDENT>else:<EOL><INDENT>yout = np.zeros( (len(t), Ny), np.float_)<EOL><DEDENT>yout[<NUM_LIT:0>] = y0<EOL>i = <NUM_LIT:0><EOL>for i in np.arange(len(t)-<NUM_LIT:1>):<EOL><INDENT>thist = t[i]<EOL>dt = t[i+<NUM_LIT:1>] - thist<EOL>dt2 = dt/<NUM_LIT><EOL>y0 = yout[i]<EOL>k1 = np.asarray(derivs(y0, thist))<EOL>k2 = np.asarray(derivs(y0 + dt2*k1, thist+dt2))<EOL>k3 = np.asarray(derivs(y0 + dt2*k2, thist+dt2))<EOL>k4 = np.asarray(derivs(y0 + dt*k3, thist+dt))<EOL>yout[i+<NUM_LIT:1>] = y0 + dt/<NUM_LIT>*(k1 + <NUM_LIT:2>*k2 + <NUM_LIT:2>*k3 + k4)<EOL><DEDENT>return yout<EOL>
Integrate 1D or ND system of ODEs using 4-th order Runge-Kutta.
This is a toy implementation which may be useful if you find
yourself stranded on a system w/o scipy.  Otherwise use
:func:`scipy.integrate`.

*y0*
    initial state vector

*t*
    sample times

*derivs*
    returns the derivative of the system and has the
    signature ``dy = derivs(yi, ti)``


Example 1 ::

    ## 2D system

    def derivs6(x,t):
        d1 =  x[0] + 2*x[1]
        d2 =  -3*x[0] + 4*x[1]
        return (d1, d2)
    dt = 0.0005
    t = arange(0.0, 2.0, dt)
    y0 = (1,2)
    yout = rk4(derivs6, y0, t)

Example 2::

    ## 1D system
    alpha = 2
    def derivs(x,t):
        return -alpha*x + exp(-t)

    y0 = 1
    yout = rk4(derivs, y0, t)


If you have access to scipy, you should probably be using the
scipy.integrate tools rather than this function.