import numpy as np                       # llama a la libreria numpy (operaciones numericas)
import matplotlib.pyplot as plt          # llama a la libreria grafica
from scipy.integrate import odeint       # llama al integrador numerico 'odeint'

# vamos a definir una funcion maricial que realiza la siguente cuenta
#
# [0    1].[  y(t)  ]
# [-w^2 0] [dy(t)/dt]
#

def sistema(Y,t,omega):
    y,ydot=Y
    return ydot,-omega**2*y

t=np.linspace(0,2*np.pi,101)             # puntos donde se evalua la funcion

y0=[1,0]                                 # posicion inicial in 1 y velocidad nula
omega=2.0                                # frecuencia angular del pendulo w^2=g/l

f=odeint(sistema,y0,t,args=(omega,))     # integrador (de Euler) 
y,ydot=f[:,0],f[:,1]

plt.figure()                                    # inicio nuevo grafico
#plt.plot(y,ydot)                               # grafico de diagrama de fases
plt.plot(t,y,t,ydot)                            # grafico y vs. t e ydot vs. t
plt.title("Solucion para omega= %5g" % omega)   # titulo 
plt.show()                                      # mostrar en pantalla