Week 3 - Solving second order ODEs

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

clear all 
close all
clc

b=0.05;      %damping coefficient
g=9.81;      %gravity in m/s^2
l=1;         %length of the pendulam string in m
m=1;         %mass of the pendulum in kg
ct=1;

%initial conditions 
theta_0=[0;3];

%time_span
t_span=linspace(0,20,500);

%solve ode
[t,result]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%plotting figure
subplot(2,1,1)
plot(t,result(:,1),'color','r')
xlabel('Time')
ylabel('Plot')
hold on
subplot(2,1,2)
plot(t,result(:,2))
xlabel('Time')
ylabel('Plot')

%animation
x0=0;
y0=0;

for i=1:length(result)
    x1=l*sin(result(i,1));
    y1=-l*cos(result(i,1));
    
    figure(2)
    plot([-1 1],[0 0],'linewidth',6,'color','b');
    axis([-2 2 -2 2]);
    line([x0 x1],[y0 y1],'linewidth',2,'color','b ');
    hold on
    plot(x0,y0,'O','markersize',3,'markerfacecolor','y');
    plot(x1,y1,'O','markersize',20,'markerfacecolor','r');
    title('Pendulum dampning');
    grid on 
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;                 %incrementation of ct for saving frame 
end

movie(M)
Videofile=VideoWriter('pendulum.avi','Uncompressed AVI');
open(Videofile)
writeVideo(Videofile,M)
close(Videofile)