import numpy as np
from matplotlib import pyplot as plt

def f(t,y):
    return -100.0*(y - np.sin(t))

def ForwardEuler(t0,y0,T,h):
    N = int((T-t0)/h)
    y = np.zeros(N)
    t = np.zeros(N)
    t[0] = t0
    y[0] = y0
    for n in range(1,N):
        y[n] = y[n-1] + h*f(t[n-1],y[n-1])
        t[n] = t[n-1] + h
    return t,y

t0 = 0
y0 = 1
T  = np.pi

h  = 1.0/52.0
t,y = ForwardEuler(t0,y0,T,h)

plt.plot(t,y)
plt.xlim([-0.1,4])
plt.title("Forward Euler, h ="+str(h))

h  = 0.01
t,y = ForwardEuler(t0,y0,T,h)

plt.plot(t,y)
plt.xlim([-0.1,4])
plt.title("Forward Euler, h ="+str(h))

def BackwardEuler(t0,y0,T,h):
    N = int((T-t0)/h)
    y = np.zeros(N)
    t = np.zeros(N)
    t[0] = t0
    y[0] = y0
    for n in range(1,N):
        t[n] = t[n-1] + h
        y[n] = (y[n-1] + 100.0*h*np.sin(t[n]))/(1.0 + 100.0*h)
    return t,y

def Trapezoidal(t0,y0,T,h):
    N = int((T-t0)/h)
    y = np.zeros(N)
    t = np.zeros(N)
    t[0] = t0
    y[0] = y0
    for n in range(1,N):
        t[n] = t[n-1] + h
        y[n] = ((1.0-50.0*h)*y[n-1] + 50.0*h*(np.sin(t[n-1])+np.sin(t[n])))/(1.0 + 50.0*h)
    return t,y

h  = 0.1
t,y = BackwardEuler(t0,y0,T,h)

plt.plot(t,y)

t,y = Trapezoidal(t0,y0,T,h)
plt.plot(t,y)

plt.xlim([-0.1,4])
plt.title("Step size h ="+str(h))
plt.legend(("Backward Euler","Trapezoid"))