def rk4_sys(f,ic,t0,T,stepsize):

    # A routine to implement RK4 method for a system (also works on a scalar equation).
    # Key points:
    # The spatial variables are y[0], y[1], ...,y[d-1], d deduced from initial data "ic".
    # t0 is the initial time, T is the final time
    # The vector "f" is a funciton involving t,y0,...,y[d-1] that defines y'=f(t,y)
    # Solution times are returned in row vector "tv"
    # Solution values are returned in matrix "yv". Rows correspond to times, columns
    # to variables y1,...,yd.

    #Spatial dimension deduced from initial data.
    d = len(ic);

    #Convert all to numeric, just to be safe
    t0 = t0.n(); T = T.n(); stepsize = stepsize.n();

    # Bookkeeping, we'll take Nf steps. Also define arrays tv and yv for time, solution.
    N = floor((T-t0)/stepsize);
    T1 = t0+N*stepsize;
    Nf = N+1;
    if abs(T1-T) < 1.0e-8:
        Nf = N;
    tv = vector(SR,Nf+1);
    yv = matrix(SR,Nf+1,d);

    #Set initial condition
    tv[0] = t0;
    for i in range(d):
        yv[0,i] = ic[i];

    vv = vector(SR,d);
    # March solution out in time.
    for i in range(1,Nf):
        told = tv[i-1];
        tv[i] = told + stepsize;
        yold = yv.row(i-1);
        m1 = f(told,yold);
        for j in range(0,d):
            vv[j] = yold[j]+0.5*stepsize*m1[j];
        m2 = f(told+0.5*stepsize,vv);
        for j in range(0,d):
            vv[j] = yold[j]+0.5*stepsize*m2[j];
        m3 = f(told+0.5*stepsize,vv);
        for j in range(0,d):
            vv[j] = yold[j]+stepsize*m3[j];
        m4 = f(told+stepsize,vv);
        m = (m1+2*m2+2*m3+m4)/6.0;
        for j in range(0,d):
            yv[i,j] = yv[i-1,j] + stepsize*m[j];

    #Last special step
    told = tv[Nf-1];
    stepsize = T - told;
    tv[Nf] = T;
    yold = yv.row(Nf-1);
    m1 = f(told,yold);
    for j in range(0,d):
        vv[j] = yold[j]+0.5*stepsize*m1[j];
    m2 = f(told+0.5*stepsize,vv);
    for j in range(0,d):
        vv[j] = yold[j]+0.5*stepsize*m2[j];
    m3 = f(told+0.5*stepsize,vv);
    for j in range(0,d):
        vv[j] = yold[j]+stepsize*m3[j];
    m4 = f(told+stepsize,vv);
    m = (m1+2*m2+2*m3+m4)/6.0;
    for j in range(0,d):
            yv[Nf,j] = yv[Nf-1,j] + stepsize*m[j];

    return tv,yv