Week 3 - Solving second order ODEs

%Solving second order ODE of simple pendulum

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

%initial condition
t0=[0;3];

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

%solving ODE
[t,results]=ode45(@(t,theta) ode_fun(t,theta,b,g,l,m),t_span,t0);

%plotting
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('Plot')
theta1=results(:,1);

%pendulum animation
c=1;
for i=1:length(theta1)
    theta01=theta1(i);
    %cordinates
    x0=0;
    y0=0;
    x1=-l*sin(theta01);
    y1=-l*cos(theta01);

    %plotting
    figure(2)
    plot([x0,x1],[y0,y1],LineWidth=2)
    %to mark sphere at the end of pendulum
    hold on
    plot(x1(end),y1(end),Marker="o",MarkerSize=20,MarkerFaceColor="r")  
    hold off
    axis([-2 2 -2 2])
    grid on
    pause(0.03)
    anim(c)=getframe(gcf);
    c=c+1;
end

%for getting animated video file
movie(anim)
animfile=VideoWriter('simple pendulum.avi','Uncompressed AVI');
open(animfile)
writeVideo(animfile,anim)
close(animfile)

%ODE function definition
function [dtheta_dt]= ode_fun(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