from __future__ import division
from pylab import *
from scipy.optimize import fsolve

def euler(f, t0, y0, h, N):
    t = t0 + arange(N+1)*h
    y = zeros((N+1, size(y0)))
    y[0] = y0
    for n in range(N):
        y[n+1] = y[n] + h*f(t[n], y[n])
    return y

def backwardEuler(f, t0, y0, h, N):
    t = t0 + arange(N+1)*h
    y = zeros((N+1, size(y0)))
    y[0] = y0
    for n in range(N):
        def F(ynplus1):
            return -ynplus1 + y[n] + h*f(t[n+1], ynplus1)
        y[n+1] = fsolve(F, y[n])
    return y