clear allclose allclcb = 0.05g = 9.81l = 1m = 1theta_0 = [0;5]t_span = linspace(0,20,500)[t,results] = ode45( @(t,theta) ode_func1(t,theta,b,g,l,m),t_span, theta_0);subplot(5,4,[13,14,17,18]);plot(t,results(:,1),'color','b');hold onplot(t,results(:,2),'color','r');hold offsubplot(5,4,[15,16,19,20]);plot(results(:,1), results(:,2));xlabel('position-Radions');ylabel('Angular velocity'); ct =1;for i =1:length(results(:,1)); x0 = 0; y0 = 0; x1 = 1*sin(results(i,1)); y1 = -1*cos(results(i,1)); subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]); plot([-1 1],[0 0],'linewidth',4,'color','k'); hold on line([x0 x1],[y0 y1],'linewidth',2,'color','b'); hold on; r=0.075; cir_deg =linspace(0,2*pi,1000); xc=r*cos(cir_deg); yc=r*sin(cir_deg); plot(x1+xc,y1+yc,'linewidth',2)plot(x1,y1)r = 0.25; rectangle('Position',[(x1-(r/2)),(y1-r/2),r,r],'Curvature',[1 1],'Facecolor','r','linewidth',3)hold onaxis equalgrid on;axis([-1.5 1.5 -1.5 0.5])pause(0.0003);hold offq(ct) = getframe(gcf);Text = sprintf("pendulum%e.png",ct);ct = ct+1;saveas(1,Text)endmovie(q)videofile = VideoWriter('Pendulum_damped.avi','uncompressed avi')open(videofile)writeVideo(videofile,q)close(videofile)Sub Function ode_func1function [dtheta_dt] = ode_func(t,theta,b,g,l,m);theta1 = theta(1);theta2 = theta(2);dtheta1_dt = theta2;dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);dtheta_dt = [dtheta1_dt; dtheta2_dt];end