solving second order ode and pendulum animation

% solvind ode and pendulum animation 
clc
%input values
b=0.5;
g=9.81;
l=1;
m=1;

%initial conditions
initial = [0;3];
% time points
t_span = linspace(0,20,500);

% solving given ode
[t,results] = ode45(@(t,theta) odefunction(t,theta,b,g,l,m),t_span,initial);
% plotting graphs
figure(1)
plot(t,results(:,1))
hold on
grid on
plot(t,results(:,2)) , xlabel(\'time\'), ylabel(\'plot\')

% creating animation of pendulum
figure(2)
grid on
axis([-0.05 0.05 -1.5 0.1]);   
origin = [0 0];
ct=1;
for i=1:length(t)
  
p = [l*sind(results(i,1)) l*-cosd(results(i,1))];
o_circ = viscircles(origin,0.001);

pendulum = line([origin(1) p(1)],[origin(2) p(2)]);
q = [p(1) p(2)];
ball = viscircles(q,0.005);
pause(0.001);
M(ct) = getframe(gcf);
ct=ct+1;

if i<length(t)
    delete(pendulum)
    delete(ball)
    delete(o_circ)
end
end

movie(M)
videofile = VideoWriter(\'pendulum4\');
open(videofile)
writeVideo(videofile,M)
close(videofile)

% function coding

function dtheta_dt = odefunction(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];

% link for pendulum animation

https://youtu.be/-ikarvZnddU