Solving second order ODE of simple pendulum and Simulating its motion

1. Link for the animation

https://youtu.be/CcEEDS8ZekU

2.Code

% Initial commands 
close all;
clc
clear all;

%--------------------------------------------------------------------------
%% Giving the dimensions ans also required parameters for the pendulum
%--------------------------------------------------------------------------
x0=0;
y0=0;
%theta= [0:10:90];
theta_I=[0;20];
len=2;
c=.01;%Damping Coefficient
g=9.81;
m=.2;
r=.2;%radius of the bob
%--------------------------------------------------------------------------
%% Defining the function for the simple pendulum and the required parameters fos it.
%--------------------------------------------------------------------------
time_span= linspace(0,150,600);%Defining the run time of the programe
[t,theta_final] = ode45(@(t,theta) P(t,theta,c,m,g,len),time_span,theta_I);
figure(1)
plot(t,theta_final(:,1))%angular accleration
hold on
plot(t,theta_final(:,2))%angular displacement

%--------------------------------------------------------------------------
%% Generation and plotting of the pendulum
%--------------------------------------------------------------------------
ct=1;
for i=1:length(theta_final(:,1));
    %The position of the bob
    x1=sind(theta_final(i,1)+180)*len;
    y1=cosd(theta_final(i,1)+180)*len;
   figure(2)
    plot([x0,x1],[y0,y1],\'linewidth\',2)
    %circles(x1,y1,\'r\',\'color\',\'g\')
    viscircles([x1,y1],r)
    axis([-2 2 -3 0])
    M(ct)=getframe(gcf);
    ct=ct+1;
end
%--------------------------------------------------------------------------
%% Creation of the Uncompressed Avi video
%--------------------------------------------------------------------------
movie(M)
% At first  create a file and specify the kind of formate required
videofile1=VideoWriter(\'Simple_pendulum.avi\',\'Uncompressed AVI\');
%Open the file that was created
open(videofile1)
%Writing all the image files into a video 
writeVideo(videofile1,M)
close(videofile1)



% The function

function [ theta_out ] = P( t,theta,c,m,g,l )
%--------------------------------------------------------------------------
%% Function describing the 2nd order differential equation for the Simple Pendulum
% C= Damping factor
% m= mass of the bob
% l= length of the pendulum
%--------------------------------------------------------------------------
k=(m*g)/l;
dam_rat=(c/(2*sqrt(m*k)));
omega=sqrt(k/m);
theta1=theta(1);
theta2=theta(2);
dtheta1_dt=theta2;
dtheta2_dt=-2*dam_rat*omega*dtheta1_dt-omega^2*sind(theta1);
theta_out=[dtheta1_dt;dtheta2_dt];

end