Solving second order ODEs in Python

CODE:

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.02
g=9.81
l=1
m=0.1

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

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

# SOLVE ODE
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\\theta_2}{dt}=-frac{b}{m}\\theta_2-\\frac{g}{l}sin\\theta1$\')
plt.ylabel(\'Plot\')
plt.xlabel(\'time\')


plt.legend(loc=\'best\')
plt.show()

#ANIMATION
ct=1
for Beta in theta[:,0]:
	x0=0
	y0=0
	x1= 1*math.sin(Beta)
	y1= -1*math.cos(Beta)
	


	filename=str(ct)+\'.png\'
	ct=ct+1

	plt.figure()
	plt.title(\'Animation of the given pendulum bob\')
	plt.plot([0,0],[1,-1],\'r--\')
	plt.plot([-1,1],[0,0],\"brown\",linewidth=10) #Create the base
	plt.plot([x0,x1],[y0,y1],\"black\") #Creates the string of the pendulum
	plt.plot(x1,y1,\"o\",markersize=20) #Creates the pendulum bob
	plt.xlim([-1.5,1.5]) 
	plt.ylim([-1.5,1.5])
	plt.savefig(filename)
	

Link to the video of animation on youtube : https://youtu.be/Kn6ojd3dUlU

Errors :

Solution :

Clear memory, close app and reopen.