Solving second order ODEs

%Calling a function - pendulum ode solver

function [dtheta_dt] = func_theta(t,theta,g,l,m)

theta1 = theta(1);
theta2 = theta(2);

dtheta1_dt=theta2;
dtheta2_dt=-((g^2)/(l^2*m^2))*sin(theta1);
dtheta_dt=[dtheta1_dt ; dtheta2_dt];
end

%Program

clear all
close all
clc

%input
g=9.81;
l=1;
m=0.6;

theta_0=[0;20];
t_span=linspace(0,10,500);

[t,result]=ode45(@(t,theta) func_theta(t,theta,g,l,m),t_span,theta_0);
figure(1)
hold on
plot(t,result(:,1),\'linewidth\',1.5\')

plot(t,result(:,2),\'linewidth\',1.5\')

figure (2)

alpha = result(:,1);

Xa1=2;
Ya1=2;

Xa2=-2;
Ya2=2;

%Creating loop

ct=1;

for i=1:length(t)

x0 = 0;
y0 = 2;

x1 = l*sin(alpha(i));
y1= -l*cos(alpha(i));

r = 0.2;
xm = x1;
ym = y1;

theta_m = linspace(0,2*pi);
x_m = r*cos(theta_m) + xm;
y_m = r*sin(theta_m) + ym;

plot([x0 x1],[y0 y1],[Xa1 Xa2],[Ya1 Ya2],x_m,y_m,\'linewidth\',3)

% plot([Xa1 Xa2],[Ya1 Ya2])

axis([-3 3 -3 3])
pause(0.03)
grid on

M(ct)= getframe(gcf);

ct=ct+1;

end

movie(M)
videofile = VideoWriter(\'ode_pen.avi\',\'Uncompressed AVI\');
open(videofile)
writeVideo(videofile,M)
close(videofile)

video link: https://youtu.be/bkma4nz3ZvA