Week 3 - Solving second order ODEs

Code for pendulum motion simulation:

%Simple pendulum

%Damping coeff
b=0.05;
%Acceleration due to gravity in m/s2
g=9.81;
%Length in m
l=1; 
%Mass in kg
m=1;

%Initial conditions - displacement & Velocity
theta_0=[0;3];

%Time span
t_span=linspace(0,20,500);

%Sloving ODE
[t, result]=ode45(@(t,theta) ode_func_sp(t,theta,b,g,l,m),t_span, theta_0);

%Pendulum motion simulation
ct=1;
for i=1:length(result(:,1))
    x0=1; y0=1;
    x1=l-sin(result(ct,1));
    y1=l-cos(result(ct,1));
    figure(1)
    plot([x0 x1],[y0 y1],'linewidth',2,'Marker','o','MarkerSize',10)
    axis([0 2 -0.5 1.5])
    title('Pendulum Motion')
    pause(0.005)
    M(ct)=getframe(gcf);
    ct=ct+1;
end

%Creating video file of pendulum motion
movie(M)
videofile = VideoWriter('pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

%Plotting displacement and velocity agaist time
figure(2)
plot(t,result(:,1))
hold on
plot(t,result(:,2))
xlabel('Time(s)')
ylabel('Displacement & Velocity')
legend('Displacement(Radians)','Velocity(rad/s)')
title('Displacement & Velocity VS Time')

ODE function definition: 

%ODE function definition
function [dtheta_dt]=ode_func_sp(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