Week 4.1 - Solving second order ODEs

clear all
close all
clc

%Given data
b = 0.05;   %damping coefficient
g = 9.81;   %Acceleration due to gravity
l = 1;      %Length of Pendulum in meters
m = 1;      %Mass 

%Time duration for simulation 
t_span = linspace(0,20,100);
%Providing Initial Condition as Angular displacement & Velocity
y_initial = [0;3];

%Calculation using ODE45
a = @(t,y)ode_second_order_pendulum(y,b,m,l,g);
[t,y]= ode45(a,t_span,y_initial);

%Plot of Angular displacement & Velocity with respect to Time function.
figure(1);
grid on
subplot(1,2,1)
plot(t,y(:,1));
title('Displacement vs Time');
xlabel('Time in seconds');
ylabel('Angular displacement in rad');
subplot(1,2,2)
plot(t,y(:,2));
title('Angular velocity Vs Time')
xlabel('Time in seconds');
ylabel('Angular velocity in rad/sec');
grid off
saveas(figure(1),'subplot','jpg')

%Polar Co-ordinate To Cartesian Co-ordinate
y1 = -1*cos(y(:,1));
x1 = 1*sin(y(:,1));

%Video for Pendulum motion
for i = 1:100
    figure(2)
    plot([0 x1(i)],[0 y1(i)],'linewidth',3,'color','r');
    hold on
    plot(x1(i),y1(i),'marker','o','linewidth',10,'markerfacecolor','g');
    plot([-1.2 1.2],[0 0],'linewidth',3,'color','b');
    axis([-1.2 1.2 -1.2 1.2]);
    title(['Time: ', num2str(t_span(i))]);
    M(i) = getframe(gcf);
    close(2)
end

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

The Function is given as:

function ydot = ode_second_order_pendulum(y,b,m,l,g)
ydot = [y(2);-(b/m)*y(2)-(g/l)*sin(y(1))];
end