def f(x,y):
    return (x-y)/2  #function dy/dx = (x-y)/2

x0 = 0.0            #initial value
y0 = 1.0            #initial value
def f(x0,y0):       #initial condition
    1.0
a = 0
b = 3.0             # find for y(3) given y(0)=1- initial condition
h = 0.5
n = (b-a)/h
k1 = []
k2 = []
k3 = []
k4 = []
x = [n]
y = [n]
for i in range(0,n+1):      #iteration steps
    k1 = h*f(x[i],y[i])
    k2 = h*f(x[i] + 0.5*h,y[i] + k1*0.5)
    k3 = h*f(x[i] + 0.5*h,y[i] + K2*0.5)
    k4 = h*f(x[i] + h,y[i] + k3)
    y.append(y + 0.16*(k1 + 2.0*k2 + 2.0*k3 + k4))  #y(i+1) = y(i) + 0.16*(k1 + 2.0*k2 + 2.0*k3 + k4)
    print "y"