Week 3 - Solving second order ODEs

close all
clear all
clc


b = 0.5;
g = 9.81;
l = 0.75;
m = 1;

% initial condition
theta_0 = [0;5];

%time points
t_span =linspace(0,10,150);
    
% solve ODE
[t, results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

plot(t,results(:,1))
hold on

plot(t,results(:,2))

xlabel('plot')
ylabel('time')

ct = 1;
for i = 1:length(t_span)
    THETA = results(i,1);
    x = 0;
    y = 0;
    x1 = l*sin(THETA);
    y1 = -l*cos(THETA);

    % plotting
    figure(2)
    plot([x, x1],[y, y1],'linewidth',3)
    hold on
    plot(x1, y1, 'o')
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.03)

    M(ct) = getframe(gcf);
    ct = ct+1;
    hold off
end
movie(M)
videofile = VideoWriter('ODE45.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