def euler_method(f,ic,stepsize,finalT):
    #Basic 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, forward or backward
    for i in range(1,Nf):
        sol[i,0] = sol[i-1,0]+stepsize
        sol[i,1] = sol[i-1,1] + stepsize*f2(sol[i-1,0],sol[i-1,1])

    #Last step, might be special.
    sol[Nf,0] = tf
    sol[Nf,1] = sol[Nf-1,1] + (tf-sol[Nf-1,0])*f2(sol[Nf-1,0],sol[Nf-1,1])

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

    return sol