pendulum ODE

link for animation: https://www.youtube.com/watch?v=TaeTlK4WMK4

clear all
close all
clc
b = 0.5;
g = 9.81;
l = 1;
m = 1;
%initial conditions
theta_0 = [0;3];
%time points
t_span = linspace(0,20,1000);
%solve ODE
[t, results] = ode45(@(t,theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);
%%animation
%origin
O = [0 0];
rectangle(\'Position\', [-0.5 0 1 0.05]);
axis([-1.4 1.4 -1.4 0.2]);
grid on
%loop for animation
for i = 1:length(t_span)
    %mass point
    P = l*[sin(results(i,1)) -cos(results(i,1))];
    %circle in origin
    O_circ = viscircles(O,0.01);
    %pendulum
    pend = line([O(1) P(1)],[O(2) P(2)],\'color\',\'m\');
    %ball
    ball = viscircles(P,0.1);
    pause(0.001);
    if i<length(t)
        delete(pend);
        delete(ball);
        delete(O_circ);
    end
end
%plots
plot(t,results(:,1));
hold on
plot(t,results(:,2));
grid on

%function
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

pendulum plot: