Week 4.1 - Solving second order ODEs

clear all
close all
clc
%Inputs
b=0.05; %Damping Coefficient
g=9.81; %Gravitiy [m/s2]
l=1;    % Length of Pendulum [m]
m=0.1;  %Mass of ball [Kg]
ct=1;

%inital condition [angular displacment, Velocity[rad/s]
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);
%plot
%color order for y axises 
colororder({'r','b'})
%left y axis
yyaxis left
%displacement plot
plot(t,results(:,1),'Color','r','LineStyle','-')
xlabel('time [s]')
ylabel('Displacment [m]','Color','r')
hold on
%velocity plot
plot(t,results(:,2),'Color','b','LineStyle','-')
yyaxis right
ylabel('Velocity [rad/s]','Color','b')
legend('Displacment [m]','Velocity [rad/s]')                               %legend
title('Pendulum ODE','FontSize',11)                         %title 

% Pendulum simulation
% Displacement results
   theta2 = results(:,1);
%for loop simulate the position of the pendulum

for i = 1:length(theta2)
    
    theta_x = theta2(i);         
    %pendulum pivot point(0,0)
    x0= 0;
    y0= 0;
    %pendulum end point effected by theta value
    x1 = l*sin(theta_x); 
    y1 = -l*cos(theta_x);
   
    
    %plotting the pendulum
    figure(2)
    plot([-0.5 0.5],[0 0],'color','g','linewidth',10);
    hold on 
    plot([x0 x1],[y0 y1],'color','r','linewidth',2)
    hold on 
    plot(x1,y1,'o','markers',15,'markerfacecolor','b')
    hold off
    grid on 
    axis([-1.5 1.5 -1.5 0.5]);
    
    pause(0.002)
    %capturing each frame 
    M(ct) = getframe(gcf);
    %changing ct index
    ct = ct+1;
end
%animation 
movie(M)
videofile = VideoWriter('ODE_pendulum.avi','Uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)


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

theta1=theta(1);
theta2=theta(2);
%first ODE
dtheta1_dt=theta2;
%Second ODE
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
%Solution Matrix
dtheta_dt=[dtheta1_dt;dtheta2_dt];

Pendulum Animation Link- https://youtu.be/ey-NCT1qQaA