Week 3 - Solving second order ODEs

clear all
close all
clc

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

%initial conditions
theta_0=[0;3];

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

[t,theta]= ode45('ode_func1', t_span, theta_0);

figure(1)
plot(t,theta(:,1))
hold on
plot(t,theta(:,2),'red')
legend('displacement vs time','velocity vs time')

figure(2)
frame=1;
for i=1:length(t_span)
    disp(i)=theta(i,1);

a(i)=l*sin(disp(i));
b(i)=-l*cos(disp(i));

x_cord(i)=a(i);
y_cord(i)=b(i);

plot([0,x_cord(i)],[0,y_cord(i)])
line([-0.5 0.5],[0 0],'linewidth',2)
ball=viscircles([x_cord(i),y_cord(i)],0.05);
grid on
axis([-2 2 -2 1])
pause(0.00001)
M(frame)=getframe(gcf);
frame=frame+1;
end

movie(M)
videofile=VideoWriter('Simple_Pendulum.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

ode_func1 code

function [dtheta_dt]= ode_func1(t,theta,b,g,l,m)
b=0.5;
g=9.81;
l=1;
m=1;
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