from pylab import *
from numpy import array,zeros,arange
def oscilator(state,time):
    k=2
    m=1
    v0=state[1]
    a=-k/m*state[0]
    return array([v0,a])
def Euler(y,t,dt,derivative):
    y_next=y+(derivative(y,t)+derivative(y+derivative(y,t)*dt,t+dt))*dt/2.
    return y_next
xo=10 # initial conditions
vo=0
dt=0.1
to=0
te=500
t=arange(to,te,dt)
N=len(t)
y=zeros([N,2])
y[0,0]=xo
y[0,1]=vo

for i in range(N-1):
    y[i+1]=Euler(y[i],t[i],dt,oscilator)

position=[y[j,0] for j in range(N)]
velocity=[y[j,1] for j in range(N)]
subplot(2,1,1)
plot(t,position)
ylabel("x")
xlabel("time")
subplot(2,1,2)
plot(t,velocity)
ylabel("v")
xlabel("time")
show()