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

## definition of the differential system
## We set the natural frequency omega0 to 1

def pendulum(X, t):
    return [X[1], - np.sin(X[0])]

omega0=1
T0= 2*np.pi/omega0 # period of linear oscillations
t = np.linspace(0, 2.5*T0, 1000) # maximum t must be sufficiently large to obtain closed orbits for all initial conditions.

###### define of set of initial conditions  in the interval (0,pi)  #####

theta0_init = np.linspace(0, np.pi, 20)

plt.figure(figsize=(6, 6))

for theta0 in theta0_init:
    Xvalues = odeint(pendulum, (theta0, 0), t)
    plt.plot(Xvalues[:,0], Xvalues[:,1])

plt.xlabel(u'$\\theta$', fontsize=26)
plt.ylabel(u'$\dot{\\theta}$', fontsize=26)
plt.axis('equal')
plt.xlim(-np.pi, np.pi)
plt.ylim(-2, 2)
plt.title('Phase Diagram')
plt.tight_layout()
plt.show()