Week 3 - Solving second order ODEs

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

#fuctions that returns dz/dt
def model(theta,t,b,g,l,m):
    
    theta1 = theta[0]
    theta2 = theta[1]
    dtheta1_dt = theta2
    dtheta2_dt = -(b/m)*theta2-(g/1)*math.sin(theta1)
    dtheta_dt = [dtheta1_dt,dtheta2_dt]    
    return dtheta_dt

b=0.05
g=9.81
l=1
m=1

#initial condition
theta_0 = [0,5]

#time points
t = np.linspace(0,10,150)

#solving ode
theta = odeint(model,theta_0,t,args = (b,g,l,m))

 #plotting
plt.plot(t,theta[:,0],'b-',label=r'$\frac{d\theta_1}{dt}=\theta2$')
plt.plot(t,theta[:,1],'r--',label=r'$\frac{d\theta2}{dt}=\frac{b}{m}\theta_2-\frac{g}{l}sin\theta_1$')
plt.ylabel('plot')
plt.xlabel('time')
plt.legend(LOC='best')
plt.show()