Week 3 - Solving second order ODEs

%ode_function

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

%inputs

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

%initial condition
theta_0 = [0;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);
plot(t, results(:,1))
hold on
plot(t, results(:,2))
ylabel('plot')
xlabel('time in sec')

% animation of simple pendulum
ct = 1;
theta_1 = results(:,1);
for i = 1:length(results(:,1))
   %initial points
   x_start = 0;
   y_start = 0;
   theta_2 = theta_1(i);
   
   %final points
   x_end = 1*sin(theta_2);
   y_end = -1*cos(theta_2);
   
   %plot for animation
   figure(2)
   %base of pendulum
   plot([-1,1],[0,0],'linewidth',4,'color','black')
   hold on
   %string of pendulum
   line([x_start x_end],[y_start y_end],'linewidth',6,'color','b')
   hold on
   %bob of pendulum
   plot(x_end,y_end,'marker','o','markersize',28,'markerfacecolor','r')
   axis([-2 2 -2 2])
   pause(0.03)
   hold off
   M(ct) = getframe(gcf);
   ct =ct+1;
end
%storage as a movie
movie(M)
videofile = VideoWriter('Ode_pendulum','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);