Week 3 - Solving second order ODEs

the program is written for the ODE of a pendulum which is attached by a string of length 'l' have mass 'm', damping coefficient 'b', and g=9.81 m/s2. . By adding forces and torques to this model, you incrementally change the pendulum from undamped and free to damped and driven. the forces acting are gravitational force, the force of inertia, and damping force, then solving the equation using the ODE function and plotting the graph of results concerning the time.

% solve the given ODE below

close all
clear all
clc

%given values
b = 0.05;  % damping coefficient.
g = 9.81;  % gravity in m/s^2,
l=1;       % length of the pendulum in m,
m=1;       % mass of the ball in kg,

%intial condition
theta_0 = [0;5];
ct=1
p0 = 0;
q0 = 0;
%time
t_span = linspace(0,20,500);

%solve ODE
[t,r] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span, theta_0)

for i = 1:length(t_span)
    ang_disp = r(:,1);   %angular displacement  
    p1 = l*sin(ang_disp(i));
    q1 = -l*cos(ang_disp(i));
    plot([p0 p1],[q0,q1],'linewidth',3,'Marker','o')
    hold on
    plot(p1,q1,'marker','o','Color','g','MarkerSize',30,'MarkerFaceColor','g')
    title({'Simple harmonic motion of pendulum'},'color','r')
    hold off
    axis([-2 2 -2 0])
    pause(0.03)
    xlabel('time')
    ylabel('angular displacement & angular velocity')
    M(ct) = getframe(gcf);
    ct =ct+1;
end

movie(M);
videofile = VideoWriter('pendulum.m','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile)



%function to sove the ode

function [dtheta_dt] = ode_function(t,theta,b,g,l,m)
p = theta(1)
q = theta(2)
dp_dt = q;
dq_dt = -(b/m)*q - (g/l*sin(p));
dtheta_dt = [dp_dt;dq_dt];
end

https://drive.google.com/file/d/1pGMtGSYjpFy3z1_CjsfEoJ9ZBCx5CO1T/view?usp=sharing