Week 4.1 - Solving second order ODEs

% Test Solving 2nd order ODE 
clear all 
close all
clc

% Inputs
% l = length of pendulum in m 
% m = mass of pendulum in kg
% b = damping coefficient (no unit)
% g = Gravity in m/s^2, Length in m, Mass in kg

l = 1;
m = 1;
b = 0.05;
g = 9.81;

% Initial Condition
% Initial displacement is 0
% Initial velocity is 3 m/s
theta_0 = [0;3];

% Time Points 
% t_span is variation in motion of pendulum 
t_span = linspace(0,20,500);

% Solve ODE
% ode45 is the program solver in matlab 
% ode_func is the separate code used to minimize the code 
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

plot(t,results(:,1),'color','g')
hold on 
plot(t,results(:,2),'color','r')
hold on 
subplot(5,4,[15,16,19,20]);
plot(results(:,1),results(:,2))
xlabel('Position (in radian)')
ylabel('Angular Velocity')
title('second order ode')

% To run loop for different values of theta 
ct = 1;
for i = 1:length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = 1*sin(results(i,1));
    y1 = -1*cos(results(i,1));
    subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);
    plot([1 1],[0 0],'linewidth',4,'color','k');
    hold on
    line([x0 x1],[y0 y1],'linewidth',2,'color','b');
    hold on
    
    % Making a marker shows a bob 
    plot(x1,y1,'-o','markers',20,'markerfacecolor','r');
    grid on;
    axis([-1.5 1.5 -1.5 0.3])
    pause(0.03);
    hold off
    M(ct) = getframe(gcf);
    ct = ct+1;
end
movie(M)
videofile = VideoWriter('test_second_order_ode.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)