Week 4.1 - Solving second order ODEs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Function code %%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%% Main Code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all 
close all
clc
%input variables
l=1; %length of pendulum
m=1;
b=0.5;
g=9.81;
t_span = linspace(0,20,100);
theta_0 = [0;3]; % angular velocity of 3
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span,theta_0);
%plot(t,results(:,1))
%hold on
%plot(t,results(:,2))
%xlabel('time')
%label('disp & vel')
disp = results(:,1);
d = disp.';
ct = 1;
x_d = [];
y_d =[];
%for loop for calculating the (x,y) cordinates
for i = d[1;,]; % for each value of claculated d
    fact = i/1; % i divided by length of the pendulum
    alpha = atand(fact); %calculating alpha 
    x_d(end+1)= cosd(alpha); % calculating x distance
    y_d(end+1)= sind(alpha); % calculating y distance
end
ct = 1;
%for loop for ploting the position of pendulum
for j = 1:length(x_d);
    plot ([0 x_d(j)],[0 y_d(j)],'-o','linewidth',3);
    axis([-1.5 1.5 -1.5 1.5]); %axis limit
    pause(0.03);
    M(ct) = getframe(gcf);
    ct=ct+1;
end
movie(M)
videofile = VideoWriter('pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Output : (animation AVI file attached)