second order ODE

Objective: import numpy as npfrom scipy.integrate import odeintimport matplotlib.pyplot as pltimport math # function that returns dz/dtdef model(theta,t,b,g,l,m): theta1=theta[0]theta2=theta[1]dtheta1_dt=theta2dtheta2_dt=-(b/m)*theta2-(g/1)*math.sin(theta1)dtheta_dt=[dtheta1_dt,dtheta2_dt]return dtheta_dt #b=0.02b=0.05g=9.81l=1m=0.1…