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/l)* 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_O = [0,3]

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

#solve ODE
theta = odeint(model,theta_O,t,args = (b,g,l,m))

ct = 1
filename = 'test%05d.png' % ct

#plot results
plt.figure()
plt.plot(t,theta [:,0])
plt.plot(t,theta [:,1])
plt.ylabel('plot')
plt.xlabel('time')
plt.savefig(filename)
plt.show ()