dynamical-lynch/scripts/program-02e.py

# Program 02e: Numerical and truncated series solution
# figure 2.6
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint


def ode(X, t):
    x = X[0]
    y = X[1]
    dxdt = y
    dydt = x-t**2*y
    return [dxdt, dydt]


X0 = [1, 0]
t = np.linspace(0, 10, 1000)
sol = odeint(ode, X0, t)
x = sol[:, 0]
y = sol[:, 1]
fig, ax = plt.subplots()
ax.plot(t, x, label='Numerical')
ax.plot(t, 1+t**2/2+t**4/24, 'r-', label='Truncated series')
plt.xlabel('t')
plt.ylabel('x')
plt.xlim(0, 4)
plt.ylim(0, 4)
ax.legend()
plt.show()