%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.integrate import odeint

def rhs(y,t):
    return [y[1], (F*l - m*l**2*y[3]**2*np.sin(y[2])+m*l*g*np.cos(y[2])*np.sin(y[2]))/(l*(M+m*np.sin(y[2])**2)),y[3], (F*np.cos(y[2]) -m*l*y[3]**2*np.sin(y[2])*np.cos(y[2])+(M+m)*g*np.sin(y[2]))/(l*(M+m*np.sin(y[2])**2))]

#parametros
M=5
m=1
l=1
g=9.8
F=0

# malha de tempo
t = np.arange(0,200,0.1)
# condição inicial
p0=[1,0.1,0,-0.1]

sol = odeint(rhs,p0,t)

theta = sol[:,2]
omega = sol[:,3]

plt.plot(t,theta)
plt.grid()

plt.plot(t,omega)
plt.grid()

plt.plot(theta,omega)
plt.grid()