Week 3 - Solving second order ODEs

clear all
close all
clc

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

% Inital condition
theta_0 = [0;3];

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

% solve ODE

[t, results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span, theta_0);
% Ploting result
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time t in sec');
ylabel('angular velocity & angular displacement');

ct = 1;
X0=0;
Y0=0;
% Creating for loop for animation
for i =1:length(t_span)

    X1=l*sin(results(i,1));
    Y1=-l*cos(results(i,1));% put - for pendulum swings towards down

    figure(2)% ploting with diffrent line with and color
    plot([-2 2],[0 0],'linewidth',6,'color','g');
    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('Damping of simple Pendulum');
    legend('Support rod','lenght of pendulum','Fixed pint','Bob of pendulum');
    grid on;
    hold off;
    P(ct)=getframe(gcf);
    ct=ct+1;
end
movie(P)
videofile = VideoWriter('ode_simple_pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,P)
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/2)*sin(theta1);
        dtheta_dt = [dtheta1_dt; dtheta2_dt]
end