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

    # A routine to implement Euler's 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];

    # March solution out in time.
    for i in range(1,Nf):
        tv[i] = tv[i-1] + stepsize;
        vv = yv.row(i-1);
        df = f(tv[i-1],vv);
        for j in range(0,d):
            yv[i,j] = yv[i-1,j] + stepsize*df[j];

    #Last special step
    stepsize = T - tv[Nf-1];
    tv[Nf] = T;
    vv = yv.row(Nf-1);
    df = f(tv[Nf-1],vv);
    for j in range(0,d):
        yv[Nf,j] = yv[Nf-1,j] + stepsize*df[j];

    #Return solution
    return tv,yv