Week 3 - Solving second order ODEs

% solving ODE of motion of a simple pendulum with damping
clear all
close all
clc

%inputs 
%length of pendulum(m)
l=1;
%damping coefficient
b=0.05;
%gravity in m/s2
g=9.81;
%mass of the ball in kg
m=1;

%initial condition
theta_0=[0;3];

%time span
t_span = linspace(0,20,400);
%solve ODE
[t,results] = ode45(@(t,theta)ode_fun(t,theta,b,g,l,m),t_span,theta_0);
figure(1)
plot(t, results(:,1));
hold on
plot(t, results(:,2));
xlabel('time');
ylabel('plot');
ct=1;
for i=1:length(t_span);
   x0=0;
   y0=0;
    x1=l*sin(results(i,1));
    y1=(-l*cos(results(i,1)));
    figure(2)
 
    plot([x0 x1],[y0 y1],'color','b','linewidth',3)
    hold on
    plot(x1,y1,'Marker','o','Markerfacecolor','r','Markersize',20)
    hold on
     plot([0 x1],[0 y1],'color','k','linewidth',3)
   
    plot(0, 0,'-o','MarkerFaceColor','k','Markersize',4)
    
    
    axis([-1.5 1.5 -1.5 0.5])
  grid on
    pause(0.03)
    M(ct)=getframe(gcf);
    ct=ct+1;
    hold off
end
movie(M)
videofile=VideoWriter('ode.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)


function [dtheta_dt] = ode_fun(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