# -*- coding: utf-8 -*- import numpy as np import scipy.integrate as scint import matplotlib.pyplot as plt a = 0 b = 10 w02 = 10 y0 = np.array([np.pi / 2, 0]) N = 1000 def f(y, x): return np.array([y[1], -w02 * np.sin(y[0])]) def energy(theta, thetadot): return 0.5 * thetadot ** 2 - w02 * np.cos(theta) x = np.linspace(a, b, N + 1) y_scipy = scint.odeint(f, y0, x) theta_scipy = y_scipy[:, 0] thetadot_scipy = y_scipy[:, 1] plt.plot(x, energy(theta_scipy, thetadot_scipy)) plt.plot(x, theta_scipy) plt.show()