Simple Pendulum - Solving second order ODEs

Simulation of Simple pendulum by solving 2^nd order differential equation:

A simple pendulum consists of a bob and arm fixed to a pivot point. Due to surrounding air acting as a damper, the below equation represents the motion of the pendulum.

where, b=damping coefficient(Ns/m)

m=mass (kg)

g=acceleration due to gravity(m/s^-2)

L=Length of arm(m)

theta=angle of revolution(rad)

t=time(s)

Output Images and animation:

Pendulum animation Video at below location:

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

MATLAB code:

%Clear all commands
clear all
close all
clc

%Input Definition
b=0.05;
g=9.81;
l=1;
m=1;

%initial condition
theta_t=[0;3];

%time points
t_span=linspace(0,20,500);

%solve ODE
[t,results]=ode45(@(t,theta) ode_func(t, theta, b, g, l, m), t_span, theta_t);

%Plot of Angular Displacement & Angular velocity vs Time
figure(1)
p1= plot(t, results(:,1),'b');
hold on
p2= plot(t,results(:,2),'r');
ylabel('Amplitude')
xlabel('Time in s')
hold off
legend([p1 p2], {'Angular Displacement in rad', 'Angular Velocity in rad/m^2'})

%For creating animation or simulation of pendulum
ang=results(:,1);
figure(2)

%Fix origin point
O=[0 0];

%Set aspect ratio of figure
axis(gca, 'equal');

%Set axis limit of plot, grid and title
axis([-1.2 1.2 -1.5 0.2])
grid on
title('Simulation of pendulum');

%Create or set up video
v=VideoWriter('Pendulum_Sim.avi');
v.FrameRate=60;
open(v);

%Iteration for time loop
for i=1:length(t_span)
    
    %Define point P end point of pendulum which changes respect to time
    P=l*[sin(ang(i,1)) -cos(ang(i,1))];
    
    %Define a small circle at origin to make it visible
    O_circ=viscircles(O,0.01);
    
    %Construct line between origin and variable point P for each frame
    pend = line([O(1) P(1)], [O(2) P(2)]);
    
    %Construct ball at end point of pendulum
    ball = viscircles(P, 0.05);
    
    %Give pause between loops
    pause(0.02)
    
    %Capture frame
    frame=getframe(gcf);
    writeVideo(v,frame);
   
    %Delete previous objects if it is not in final loop
    if i<length(t)
        delete(pend);
        delete(ball);
        delete(O_circ);
    end
      
end

close(v);




%ODE Function script
function[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