Week 3 - Solving second order ODEs

https://www.youtube.com/watch?v=jInBbeTtqG4

clear all
close all
clc
b=0.05;
g=9.81;
l=1;
m=1;

%Initial Conditions
theta_0=[0;3];

%Time Range
t_span=linspace(0,20,100);

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

plot(t,result(:,1),'-r')
hold on
plot(t,result(:,2),'-.b')
ylabel('plot')
xlabel('Time')
legend('Angular Displacement','Velocity')
grid on;
hold off;
%ANIMATION
X0=0;
Y0=0;
ct=1;
for i=1:length(t_span)
    X1=l*sin(result(i,1));
    Y1=-l*cos(result(i,1));
    figure(2)
    plot([-0.5 0.5],[0 0],'linewidth',7,'color','b');
    axis([-1.5 1.5 -1.5 1.5]);
    line([X0 X1],[Y0 Y1],'linewidth',3,'color','r');
    hold on;
    plot(X0,Y0,'O','markersize',3,'color','k')
    plot(X1,Y1,'.','markersize',100,'color','y')
    legend('Support Rod','Pendulum length','Fixed Point','Pendulum Bob');
    grid on;
    hold off;
    M(ct)=getframe(gcf);
    ct=ct+1;
end
movie(M);
videofile=VideoWriter('Simple Pendulum.avi');
open(videofile);
writeAnimation(videofile,M);
close(videofile);

%% Function to be Called

function [thetaprime]=f(t,theta,b,g,l,m)
theta1=theta(1);
theta2=theta(2);
thetaprime1=theta2;
thetaprime2=-(b/m)*theta2-(g/l)*sin(theta1);
thetaprime=[thetaprime1;thetaprime2]
end