from numpy import *
#import matplotlib.pyplot as plt

def euler(f, x0, t0, h, N):
    t = zeros(N+1)
    x = zeros(N+1)
    t[0] = t0
    x[0] = x0
    for k in range(N):
        x[k+1] = x[k] + h*f(t[k], x[k])
        t[k+1] = t[k] + h
    
    return x, t

def euler_mid(f, x0, t0, h, N):
    t = zeros(N+1)
    x = zeros(N+1)
    t[0] = t0
    x[0] = x0
    for k in range(N):
        x_mid = x[k] + (h/2.0)*f(t[k], x[k])
        x[k+1] = x[k] + h*f(t[k]+h/2.0, x_mid)
        t[k+1] = t[k] + h
    
    return x, t



exact = lambda t: 2*exp(t) - t - 1
f = lambda t, x: x + t


T = 1.0;
x0 = 1.0;
t0 = 0.0;

for i in range(1,6):
    N = 10**(i)
    h = T/N
    x_euler, t1     = euler(f, x0, t0, h, N)
    x_euler_mid, t2 = euler_mid(f, x0, t0, h, N)
    
    abserr_euler     = abs(x_euler[-1] - exact(T))
    abserr_euler_mid = abs(x_euler_mid[-1] - exact(T))
    print "h: %6e: err euler: %12e, err euler mid: %12e" % \
          (h, abserr_euler, abserr_euler_mid)