Week 3 - Solving second order ODEs

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

#function 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(theta_1)
	dtheta2_dt=[dtheta1_dt, dtheta2_dt]
	return dtheta_dt

# initial condition
theta_0=[math.pi/2,0]

# time points
t=np.linspace(0,15,130)

L=1
b=0.05
g=9.81
m=1
# solve
theta=odeint(model,theta_0,t,args=(b,g,L,m))

# plot results
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()