def crank_n(u0 =1.0,t0 = 0.0, T=1.0,N=100):
    dt = (1.0*(T-t0))/N
    u = [u0,]
    t = [t0,]
    
    tmp = (1. +dt*0.5)/(1. - dt*0.5)
    for i in range(N):
        u +=[u[-1]*tmp]
        t +=[t[-1]+dt]
    return u,t


def feuler(u0 =1.0,t0 = 0.0, T=1.0,N=100):
    dt = (1.0*(T-t0))/N
    u = [u0,]
    t = [t0,]
    
    tmp = (1. +dt)
    for i in range(N):
        u +=[u[-1]*tmp]
        t +=[t[-1]+dt]

    return u,t


def beuler(u0 =1.0,t0 = 0.0, T=1.0,N=100):
    dt = (1.0*(T-t0))/N
    u = [u0,]
    t = [t0,]
    
    tmp = 1.0/(1. -dt)
    for i in range(N):
        u +=[u[-1]*tmp]
        t +=[t[-1]+dt]

    return u,t

def heun(f,u0 =1.0,t0 = 0.0, T=1.0,N=100):
    dt = (1.0*(T-t0))/N
    u = [u0,]
    t = [t0,]
    
   
    for i in range(N):
        F1 = f(t[-1],u[-1])
        F2 = f(t[-1] + dt,u[-1] + dt*F1)
        u +=[u[-1] + dt*0.5*(F1 + F2)]
        t +=[t[-1]+dt]

    return u,t

def rk(f,u0 =1.0,t0 = 0.0, T=1.0,N=100):
    dt = (1.0*(T-t0))/N
    u = [u0,]
    t = [t0,]
    
   
    for i in range(N):
        F1 = f(t[-1],u[-1])
        F2 = f(t[-1] + 0.5*dt,u[-1] + 0.5*dt*F1)
        F3 = f(t[-1] + 0.5*dt, u[-1] + 0.5*dt*F2)
        F4 = f(t[-1] + dt, u[-1] + dt*F3)
        u +=[u[-1] + (dt/6.0)*(F1 + 2.0*F2 +2.0*F3 + F4)]
        t +=[t[-1]+dt]

    return u,t

from pylab import plot,show,ylabel,xlabel
from math import exp, fabs
#test
s = crank_n()
#plot(s[1],s[0])
s = feuler()
#plot(s[1],s[0])
s = beuler()
#plot(s[1],s[0])
s = heun(lambda t,x: x)
#plot(s[1],s[0])
s = rk(lambda t,x: x)
#plot(s[1],s[0])
#show()


N = [100,200,400,800,1600,3200,6400]
Tt = 5.0
f = lambda t,x: x
analytic = exp(Tt)
print "dt    C-N              FE              BE             HE              RK"
for n in N:
    dt = Tt/n
    cr = crank_n(T=Tt,N = n)
    fe = feuler(T=Tt,N = n)
    be = beuler(T=Tt,N = n)
    he = heun(f,T =Tt, N = n)
    r = rk(f,T =Tt, N = n)
    print dt," ",fabs(analytic - cr[0][-1])/(dt*dt)," ",fabs(analytic - fe[0][-1])/(dt)," ",fabs(analytic - be[0][-1])/(dt)," ",fabs(analytic - he[0][-1])/(dt*dt)," ",fabs(analytic - r[0][-1])/(dt*dt*dt*dt)