def imp_euler(f,ic,stepsize,finalT):
   #Improved Euler's method implementation. If number of steps from initial to final time is not an integer, last step is incremental of proper length.
   #Also, can solve backwards in time. Uses transformation t->-t, so u'(t) = f(t,u(t)) with u(a) = b becomes v'(t) = -f(t,v(t)) with v(-a) = b under time reversal.

    if finalT > ic[0]: #Solve forward in time
        F = 1
        t0 = ic[0]
        u0 = ic[1]
        tf = finalT
        f2 = lambda t,u: f(t,u)
    else:              #Solving backward in time
        F = -1
        t0 = -ic[0]
        u0 = ic[1]
        tf = -finalT
        f2 = lambda t,u: -f(-t,u)

    #Make sure everyone is a float
    t0 = t0.n(); u0 = u0.n(); tf = tf.n();

    N = floor((tf-t0)/stepsize) #Round down
    T1 = t0+N*stepsize #Corresponding time
    Nf = N+1 #Number of steps we'll take, unless T1 = finalT
    if abs(T1-tf) < 1.0e-8: #To some tolerance
        Nf = N
    sol = matrix(SR,Nf+1,2) #Solution time and values
    sol[0,0] = t0 #Initial time
    sol[0,1] = u0 #Initiai value

    #March out in time
    for i in range(1,Nf):
        sol[i,0] = sol[i-1,0]+stepsize #Next time
        m0 = f2(sol[i-1,0],sol[i-1,1]) #Current slope
        w = sol[i-1,1] + stepsize*m0 #A single Euler step
        m1 = f2(sol[i,0],w) #Slope at next point
        m = 0.5*(m0+m1) #Average slope
        sol[i,1] = sol[i-1,1] + stepsize*m

    #Last step, might be special.
    sol[Nf,0] = tf
    step2 = tf-sol[Nf-1,0]
    m0 = f2(sol[Nf-1,0],sol[Nf-1,1]) #Current slope
    w = sol[Nf-1,1] + step2*m0 #A single Euler step
    m1 = f2(sol[Nf,0],w) #Slope at next point
    m = 0.5*(m0+m1) #Average slope
    sol[Nf,1] = sol[Nf-1,1] + step2*m

    #If necessary, negate t values, return
    if F == -1:
        for i in range(0,Nf+1):
            sol[i,0] = -sol[i,0]

    return sol