import numpy as np

def f(Y,t):
    return np.array([Y[1] , K*w0**2-w0**2*Y[0]-2*m*w0*Y[1]])

Y0=np.array([0,0])
SY=[Y0]
h=0.001
Tmax=10
N=int(Tmax/h) #int
t=np.linspace(0,Tmax,N)
K=1
m=0.5
w0=5

#recursivité d'Euler
for i in range(0,N-1):
    SY.append(SY[i]+h*f(SY[i],t[i]))
    
SY=np.array(SY)

import matplotlib.pyplot as pl
pl.plot(t,SY[:,0])
pl.show()

import scipy.integrate as sc

SYsc=sc.odeint(f,Y0 , t)
pl.plot(t,SYsc[:,0],'r')
pl.show()