from numpy import array,arange,asarray
from pylab import plot,xlabel,ylabel,show

k=1.0 # set the spring parameter
lamb=10.0 # anharmonic parameter

def f(Y,t):
        x = Y[0]
        vx = Y[1]
        return array([vx,-k*x-lamb*x**3])

ti = 0.0
tf = 27.0
N=5000
h=(tf-ti)/N

tpoints = arange(ti,tf,h)
xpoints = []
vxpoints = []

Y = array([0.25,-0.1],float)      
for t in tpoints:
    xpoints.append(Y[0])
    vxpoints.append(Y[1])
#   This is RK4
    k1 = h*f(Y,t)
    k2 = h*f(Y+0.5*k1,t+0.5*h)
    k3 = h*f(Y+0.5*k2,t+0.5*h)
    k4 = h*f(Y+k3,t+h)
    Y += (k1+2*k2+2*k3+k4)/6
#    This is Euler
#    Y += h*f(Y,t)

plot(tpoints,xpoints)
xlabel("t")
ylabel("x(t)")
show()