Solving ODE with simple pendulum as example

This code, is to slove ODE with example of simple pendulum

clear all
close all
clc

%inputs
m=input(\'Mass of the body/shpere (m,kg):\');
l=input(\'Lenght of the string (l,m):\');
g=input(\'acceleration due to gravity (g,m/s^2):\');
b=input(\'damping coefficient (b):\');
Ad=input(\'initial angular displacement (Ad,m):\');
Av=input(\'initial angular velocity (Av,m):\');
disp(\'Enter time limits\');
Ti=input(\'initial limit (Ti,sec)\');
Tf=input(\'finial limit (Tf,sec)\');
In=input(\'no.of intervals\');
thetai=[Ad;Av];
t_limits=linspace(Ti,Tf,In);

%calling function code to solve ODE
[t,d_dt]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_limits,thetai);
figure(1);
plot(t,d_dt);
xlabel(\'Time(sec)\');
legend(\'angular displacement\',\'angular velocity\');

count=1;
for i= 1:length(d_dt(:,1))
    angulardisplacement =d_dt(i,1);
    angulardisplacement_x=+l*sin(angulardisplacement);
    angulardisplacement_y=-l*cos(angulardisplacement);
    
    %plotting
    figure(2)
    plot([-0.4 0.4],[0 0],\'k\',[-0.4 0.4],[-0.015 -0.015],\'k\',[-0.4 0.4],[0.1 0.1],\'k\',[-0.4 0.4],[0.2 0.2],\'k\');
    hold on
    plot([0 0],[0 0.2],\'k\',[-0.2 -0.2],[0 0.2],\'k\',[0.2 0.2],[0 0.2],\'k\',[0.4 0.4],[-.015 0.2],\'k\',[-0.4 -0.4],[-.015 0.2],\'k\');
    plot([0 0],[-2 0],\':k\');
    s=0.01;
    plot([(0-s) (angulardisplacement_x-s)],[(0-s) (angulardisplacement_y-s)],\'k\',[(0+s) (angulardisplacement_x+s)],[0 angulardisplacement_y],\'k\',\'linewidth\',1);
    plot(angulardisplacement_x,angulardisplacement_y,\'marker\',\'o\',\'markersize\',20,\'markeredgecolor\',\'k\',\'markerfacecolor\',\'m\');
    plot(angulardisplacement_x,angulardisplacement_y,\'.b\');
    axis([-1 1 -1.5 1]);
    hold off
    A(count)= getframe(gcf);
count=count+1;
end
movie(A)
videofile = VideoWriter(\'Simple_Pendulum.avi\',\'Uncompressed AVI\');
videofile.FrameRate = 8;
open(videofile);
writeVideo(videofile,A);
close(videofile);

ode_function code

function [dtheta_dt]= ode_func(t,theta_1,b,g,l,m)
theta1=theta_1(1);
theta2=theta_1(2);
dtheta1_dt=theta2;
dtheta2_dt= -(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

plot of angular displacemen,angular velocity w.r.t time

animation of simple pendulum