Week 3 - Solving second order ODEs

%Challenge submision -Solving second order ODEs of simple pendulum
close all
clear all
clc
% Data for pendulum
m=1.0                %mass of bob
b=0.05              %dampping coefficient
l=1.0                %lenght of string
g=9.81               %accelaration due to gravity

%Initial condition
theta_0=[0;3];%we consider initial disp. is 0 rad and initial velocity 3radian per sec   

t_span=linspace(0,20,100); %time period of 0 to 10sec with 100 points

%ode45 solver
[t_solv,theta_solv]=ode45(@(t,theta) func_dt(t,theta,b,g,l,m),t_span,theta_0);

%ploting
figure(1)
subplot(2,1,1)
plot(t_solv,theta_solv(:,1),'r','linewidth',3)
title('Angular Displacement VS Time')
xlabel('Time(sec)')
ylabel('Displacement(Radian)')
grid on;
subplot(2,1,2)

plot(t_solv,theta_solv(:,2),'b','linewidth',3)
title('Angular Velocity VS Time')
xlabel('Time(sec)')
ylabel('Angular velocity(rad/sec)')
grid on;


%Animation of Pendulum motion

x0=0;                             %initial point on support
y0=0;

theta=theta_solv                  %only for use
ct=1;                             %counert
%We use 'for loop' to creat video file and animation
for i=l:length(t_span)             
    x1=l*sin(theta(i,1));   %position of the pendulum for t1 and theta1
    y1=-l*cos(theta(i,1));
figure(2)
    plot([-1 1],[0 0],'c','linewidth',6)
    axis([-2 2 -2 2]);
    line([x0 x1],[y0 y1],'color','b','linewidth',2)
    hold on
    plot(x0,y0,'o','markersize',5,'markerfacecolor','k')
    plot(x1,y1,'o','markersize',30,'markerfacecolor','r')
    title('oscillations of simple pendulum')
    legend('Rigid support','Rope','Fixed point','Bob')
    grid on 
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
end
movie(M);
videofile=VideoWriter('Pendulum_Motion.avi')
open(videofile);
writeVideo(videofile,M);
close(videofile);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Creating function 
function dtheta_dt=func_dt(t,theta,b,g,l,m)

theta1=theta(1);         %displacement
theta2=theta(2);         %angular velocity
dtheta1_dt=theta2;
dtheta2_dt=(-b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt]
end