def implicit_euler_sys(f,Jf,ic,t0,T,stepsize):

    # A routine to implement implicit/backward 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)
    # The matrix "Jf" is the Jacobian of f with respect to 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.
    w = vector(SR,d);
    for i in range(1,Nf):
        told = tv[i-1];
        tv[i] = told + stepsize;
        yold = yv.row(i-1);
        # Need to solve w = yold + stepsize*f(told,w) for w.
        # Then Jf2(w) = I - Jf(told,w). Use Newton's method, initial guess w = yold.
        for j in range(d): #Initial guess
            w[j] = yold[j];
        for j in range(5): #5 iterations max
            r = w-stepsize*f(told,w)-yold;
            #Newton update
            ynew = w - (matrix.identity(d)-stepsize*Jf(told,w)).solve_right(r);
            if norm(ynew-w) < 1.0e-8:
                break;
            w = ynew;

        for j in range(d):
            yv[i,j] = ynew[j];

    #Here goes incremental last step.
    told = tv[Nf-1];
    tv[Nf] = T;
    stepsize = T-told;
    yold = yv.row(Nf-1);
    # Need to solve w = yold + stepsize*f(told,w) for w.
    # Then Jf2(w) = I - Jf(told,w). Use Newton's method, initial guess w = yold.
    for j in range(d): #Initial guess
        w[j] = yold[j];
    for j in range(5): #5 iterations max
        r = w-stepsize*f(told,w)-yold;
        #Newton update
        ynew = w - (matrix.identity(d)-stepsize*Jf(told,w)).solve_right(r);
        if norm(ynew-w) < 1.0e-8:
            break;
        w = ynew;

    for j in range(d):
        yv[Nf,j] = ynew[j];

    #Return solution
    return tv,yv