Week 3 - Solving second order ODEs

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

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

% initial condition
 theta_0 = [0,5];

 % TIME POINTS
  t_span = linspace(0,20,500);

 % SOLVE ODE
  [t,results] = ode45(@(t,theta) ode_func1(t,theta,b,g,l,m),t_span,theta_0);
   
   plot(t,results(:,1),'color','r')
   hold on 
   plot(t,results(:,2),'color','b')
   ylabel('displacement & angular velocity')
   xlabel('time')

% ANIMATION
ct = 1;
for i = 1:length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,l));
    y1 = -(l*cos(results(i,l)));
     
   figure(2)
   plot([-1,1],[0,0],'LineWidth',3,'color','m');
   hold on
   plot([x0 x1], [y0 y1],'LineWidth',3,'Color','b');
   hold off
   hold on
   plot(x1,y1,'o','markersize',20,'markerfacecolor','g','MarkerEdgeColor','b');
   hold off
   title('simulation of a damped pendulum motion');
   grid on
   axis([-1 1.5 -1.5 .5])

   M(ct) = getframe(gcf) ;
  ct = ct + 1;
end

% video of the animation
 movie(M)
 video=VideoWriter('PENDULUM.AVI','UNCOMPRESSED AVI');
 open(video)
 writeVideo(video,M)
 close(video)