Week 3 - Solving second order ODEs

%%%%%%%%%%%%%%%%%%%%
%%%%% ode_main.m
%%%%%%%%%%%%%%%%%%%


clear all
close all

clc

opengl('software');
opengl('save','software');

b = 0.05;
g = 9.81;
l = 1;
m = 1;
 
%initial condition
theta_0 = [0;3];

%time points
t_span = linspace(0,20,400);

%solve ode
[t, results] = ode45(@(t,theta) ode_func(t, theta, b ,g , l, m),t_span, theta_0);
figure(1);
plot(t, results(:,1),'color','g');
hold on;
plot(t,results(:,2),'color','b');
legend("angular displacement(rad)", "angular velocity(rad/s)");
ylabel("theta");
xlabel("time");
title("Motion of a simple Pendulum");
hold off
ct = 1;
for i = 1:length(t)
        x0 = 0;
        y0 = 0;
        
        x1 = l*sin(results(i,1));
        y1 = -l*cos(results(i,1));

        %plotting
        figure(2); 
        plot([x0 x1],[y0 y1],'linewidth',3);
        hold on
        plot(x1,y1,'o','MarkerSize',10,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6]);
        hold on
        plot ([-1 1],[0 0],'linewidth',5,'color','r');
        hold on
        axis([-1.1 1.1 -1.1 0])
        hold off
        M(ct) = getframe(gcf);
        ct = ct + 1;
end

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% ode_func.m
%%%%%%%%%%%%%%%%

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