def runge(de, x, tstart, tend, steps):
    t=tstart
    dt=float(tend-tstart)/steps
    for i in range(steps):
# evaluate derivative at initial point to get k1
       k1=de(t,x)
# half step forward with k1 and evaluate deriv to get k2
       k2=de(t+dt/2,[a+b*dt/2.0 for a,b in zip(x,k1)])
# half step forward with k2 and evaluate deriv to get k3
       k3=de(t+dt/2,[a+b*dt/2.0 for a,b in zip(x,k2)])
# full step forward with k3 and evaluate deriv to get k4
       k4=de(t+dt,[a+b*dt for a,b in zip(x,k3)])
# full step forward with average of k1,k2,k3,k4 for new point
       x=[a+(b1+2*b2+2*b3+b4)/6.0*dt for (a,b1,b2,b3,b4) in zip(x,k1,k2,k3,k4)]
       t+=dt
    return x

def ball(t,x):
  g=9.8
  v=x[0]
  return [g,v]



p=[0,0]
print runge(ball,p,0.0,5.0,10)
p=[0,0]
print runge(ball,p,0.0,5.0,1000)

import math
def pendulum(t,x):
  g = 9.8 #gravatational constant
  L = 1.0 #Length of pendulum

  theta = x[1]
  omega  = x[0]

  return [ -g*math.sin(theta)/L, omega]


p=[0.1,0]
dt=0.01
for i in range(1000):
  p=runge(pendulum,p,i*dt,(i+1)*dt,10)
  print p