Week 3 - Solving second order ODEs

clear all 
close all
clc

% Second order ODE

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

% initial condition
theta_0 = [0;3];
t_span = linspace(0,20,500);

% Solving ODE

[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
grid on;
axis([0 20 -3 3])


ct = 1;
theta = results(:,1);
for i = 1:length(theta)
    
% axis coordinates
x0 = 1;
y0 = 1;
x1= 1- l*sin(theta(i));
y1=-1- l*cos(theta(i));

% Plotting
figure(2)
plot([0 2], [1,1])
hold on
plot([x0 x1],[y0 y1],[x1 x1],[y1 y1],'o','MarkerSize',20,'markerfacecolor', 'r')
hold off
axis([0 2 -5 1.5])
pause(0.03);
M(ct) = getframe(gcf);
ct = ct+1;
end

movie(M)
videofile = VideoWriter('Simple_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);
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end