Week 3 - Solving second order ODEs

https://1drv.ms/w/s!AtEVXMohLim7gqRszHDmpBxMjBoCWw?e=6MBqwS

ODE_main.m file

% ODE for pendulum 

clear 
close all
clc

b = 0.05;  % damping coeff
g = 9.81;
L = 1; 
m = 1;

% I.C
theta_0 = [0;3];  % initial ang. disp. is 0 and initial ang. velocity is 3 m/s

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

% solve ODE
[t,results] = ode45(@(t,theta) ode_func(t, theta, b,g,L,m),t_span,theta_0);  
% [TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system of differential equations y' = f(t,y) 
% from time T0 to TFINAL with initial conditions Y0.




subplot(2,2,1) 
plot(t,results(:,1))  % results is disp. and velocity respectively
ylabel('Ang. Disp. (rad)')
xlabel('time (sec)')
hold on
subplot(2,2,3) 
plot(t,results(:,2))
ylabel('Ang. Velocity (rad/sec)')
xlabel('time (sec)')
hold off




% Creating an animation 

ct = 1;
for i = 1:length(results(:,1))
    % Creating the "string" for the pendulum
    x0 = 0;
    y0 = 0;
    x1 = L*sin(results(i,1));
    y1 = -L*cos(results(i,1));
    
    subplot(2,2,[2,4])
    plot([-1,1],[0,0],'color','k')
    
    line([x0,x1],[y0,y1])
    
    % Creating the filled bobber 
    radius = 0.08;
    left_x = x1 - radius;  % in (m)
    top_y = y1 - radius;
    
    % Create a circle in Matlab using the rectangle function
    rectangle('Position', [left_x, top_y, radius*2, radius*2], 'Curvature', [1,1], 'FaceColor', 'r')
    axis equal
    grid on
    axis([-1.5, 1.5, -1.5, 1.5]) 
    pause(0.1) % Frame changes every 0.1 seconds
    Figures(ct) = getframe(gcf); % "gets current figure" and saves it in an array called Figures
    ct = ct+1; % ct represents each frame

end



movie(Figures) % Matlab generates a movie of the gathered figures
videofile = VideoWriter('Pendulum_Simulation.avi','Uncompressed AVI'); % Create a new video fil(avi)
open(videofile)
writeVideo(videofile,Figures) % Write the Figures array into the file
close(videofile)

ode_func.m file (function)

function [dtheta_dt] = ode_func(t,theta,b,g,L,m)
    % theta is a variable, now we set its array values individual names
    
    theta1 = theta(1);  % Displacement
    theta2 = theta(2);  % Velocity
    
    dtheta1_dt = theta2; 
    dtheta2_dt = -(b/m)*theta2-(g/L)*sin(theta1);
    dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

https://youtu.be/s9--jPqHxqs

(Youtube vid link in word document too)