Week 3 - Solving second order ODEs

%function program
function dtheta_dt = ode_function (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

%main program
clear all
close all
clc

%given input
%damping coefficient
b = 0.05;

%gravity m/s^2
g = 9.81;

%length of pendulm in meter
l = 1;

%mass of ball in kg
m = 1

%initial condition
%angular velocity
theta_0 = [0;3];

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

%solve ode
[t, results] = ode45(@(t, theta) ode_function(t, theta,b,g,l,m), t_span, theta_0);
%plotting
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')

% Animation
x0 = 0;
y0 = 0;

ct=1;
for  i = 1 : length(t_span);
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    
    figure(2)
    
    plot([x0 x1],[y0 y1],'linewidth',3);
    axis([-2 2 -2 2]);
    hold on
    plot(x1,y1,'o','markersize',30);
    title ('simple pendulum motion');
    pause(0.03)
    hold off
    M(ct) = getframe(gcf);
    ct = ct+1;
end
%to creat movie
movie(M)
videofile = VideoWriter('simple_pendulum_motion.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M)
close(videofile);

https://youtu.be/KUcnOyV7GgE