Solving second order ODEs

%Program to solve ODE with Animation
clear all
close all
clc

%Inputs
b=0.05;
g=9.81;
l=1;
m=1;

%Initial Condition
theta_0=[0;3];

%Time Span
t_span=linspace(0,20,500);

%Solution of ODE
[t, results]=ode45(@(t,theta) pendulum(theta,b,g,l,m),t_span,theta_0);
displacement=results(:,1);

%Initial Parameters
x0=0;
y0=0;
ct=1;

%Plot
figure(1)
for i=1:length(results)
    THETA=displacement(i);
    plot([0 l*sin(THETA)],[0 -l*cos(THETA)],l*sin(THETA),-l*cos(THETA),'O')
    xlabel('X-Axis')
    ylabel('Y-Axis')
    title('Pendulum Animation')
    axis([-(l+0.2) l+0.2 -(l+0.2) 0.2])
    M(ct)=getframe(gcf);
    ct=ct+1;
end
movie(M)
videofile=VideoWriter('Pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)
figure(2)
hold on
plot(t,results(:,1))
plot(t,results(:,2))
xlabel('Time')
ylabel('Results')
title('Pendulum')
legend('Displacement','Velocity')

%Function to calculate ODE
function [dtheta_dt]=pendulum(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

https://drive.google.com/file/d/1SNrcgmcTAMV_cqWePslpQT-sjN19Fw73/view?usp=sharing