Week 3 - Solving second order ODEs

clear all
close all
clc

b = 0.05;
g=9.81;
l=1;
m=1;

%initial condition
theta_0 = [0;5];
Angular_displacement = 0
angular_velocity = 3

% time points
t_span = linspace(0,20,500);

%solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('angular_displacement vs angular_velocity')

legend ('Angular_displacement','Angular_velocity')
grid on

figure(2)
plot(results(:,1), results(:,2))
xlabel ('position')
ylabel ('Angular_velocity')
legend('Angular_velocity')

theta2 = results(:,1);
ct = 1;
 for i = 1:length(theta2)
     theta3 = theta2(i);
     % for fixed bar and movung bar of a pendulum
     x0 = 0;
     y0 = 0;
     
     x1 = 3*sin(theta3);
     y1 = -3*cos(theta3);
     
     figure(3)
     plot([-3.3 3.3],[0,0],'linewidth',4)
     hold on
     plot([x0 x1],[y0 y1], 'g','linewidth',4)
     hold on
     plot(x1,y1,'o','markers',30,'markerfacecolor','m')
     xlabel('x-axis')
     ylabel('y-axis')
     hold off
     grid on
     axis ([-3.5 3.5 -3.5 0.5])
     % time frame to get the the video
     pause(0.003)
     M(ct) = getframe(gcf);
     ct = ct+1;
 end
 
 movie(M)
 videofile = VideoWriter('Simple_Pendulum.avi','unompressed AVI');
 open(videofile)
 writeVideo(videofile,M)
 close(videofile)
 


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