import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def f(x,y):
    return 10*np.exp(-(x-2)**2/(2*(0.075**2))) - 0.6*y

h = 0.1
x0 = 0
y0 = 0.5
xlist = []
xlist.append(x0)
ylist = []
ylist.append(y0)
xi = x0
yi = y0
for i in xrange(40):
    k1 = f(xi,yi)
    k2 = f(xi+h/2,yi+k1*h/2)
    k3 = f(xi+h/2,yi+k2*h/2)
    k4 = f(xi+h,yi+k3*h)
    xi += h
    yi += (k1+2*k2+2*k3+k4)*h/6
    xlist.append(xi)
    ylist.append(yi)

h = 0.01
x0 = 0
y0 = 0.5
xlist = []
xlist.append(x0)
ylist = []
ylist.append(y0)
xi = x0
yi = y0
for i in xrange(400):
    k1 = f(xi,yi)
    k2 = f(xi+h/2,yi+k1*h/2)
    k3 = f(xi+h/2,yi+k2*h/2)
    k4 = f(xi+h,yi+k3*h)
    xi += h
    yi += (k1+2*k2+2*k3+k4)*h/6
    xlist.append(xi)
    ylist.append(yi)

h = 0.001
x0 = 0
y0 = 0.5
rk4xlist = []
rk4xlist.append(x0)
rk4ylist = []
rk4ylist.append(y0)
xi = x0
yi = y0
for i in xrange(4000):
    k1 = f(xi,yi)
    k2 = f(xi+h/2,yi+k1*h/2)
    k3 = f(xi+h/2,yi+k2*h/2)
    k4 = f(xi+h,yi+k3*h)
    xi += h
    yi += (k1+2*k2+2*k3+k4)*h/6
    rk4xlist.append(xi)
    rk4ylist.append(yi)

h = 0.001
x0 = 0
y0 = 0.5
rk5xlist = []
rk5xlist.append(x0)
rk5ylist = []
rk5ylist.append(y0)
xi = x0
yi = y0
for i in xrange(4000):
    k1 = f(xi,yi)
    k2 = f(xi+h/4,yi+k1*h/4)
    k3 = f(xi+h/4,yi+k1*h/8+k2*h/8)
    k4 = f(xi+h/2,yi-k2*h/2+k3*h)
    k5 = f(xi+3*h/4,yi+3*k1*h/16+9*k4*h/16)
    k6 = f(xi+h,yi+(-3*k1*h+2*k2*h+12*k3*h-12*k4*h+8*k5*h)/7)
    xi += h
    yi += (7*k1+32*k3+12*k4+32*k5+7*k6)*h/90
    rk5xlist.append(xi)
    rk5ylist.append(yi)

xtemp = np.asarray(rk5xlist)-np.asarray(rk4xlist)
ytemp = np.asarray(rk5ylist)-np.asarray(rk4ylist)
plt.plot(np.asarray(rk5xlist),ytemp)
# difference in the estimation of the function value using rk4 and rk5 methods with same step size, 0.001