close all
clear all
clc

vinf=5;
g=9.81;
v0=50;
f=@(v,a) g*(1-(v/vinf).^a);

alpha=1.1;

for j=1:4 % cycle to repeat the experiment for the 4 values of alpha

dt=0.05;

for i=1:3 % cycle to repeat the experiment for the 3 values of h: 0.05,0.025,0.0125

v=e73(dt,f,v0,alpha); % call to explicit euler function, specifying alpha 

% % uncomment this if you want to plot the solution
% t=0:dt:1;
% figure
% plot(t,v)
% xlabel('t')
% ylabel('v(t)')
% grid on
% title(['Speed of the parachutist as a function of time, dt=' num2str(dt) ', alpha=' num2str(alpha) ])

disp(' ')
disp(['dt=' num2str(dt) ', alpha=' num2str(alpha)])
disp(['Finale value of v: ' num2str(v(end)) ])

dt=dt/2;% halving the time-step
end
disp(' ')

alpha=alpha+0.2;% increasing alpha

end

%%

Vend1=[8.6894, 8.9888, 9.1377]; % solutions from the previous part
Vend2=[6.4882,6.7136,6.8284];
Vend3=[5.5998,5.7491,5.8273];
Vend4=[5.2417,5.3386,5.3899];

% extrapolation to get second order error from first order approximations
v1=[2*Vend1(2)-Vend1(1), 2*Vend1(3)-Vend1(2)] % case alpha=1.1

abs(v1(1)-v1(2))

%%
v2=[2*Vend2(2)-Vend2(1), 2*Vend2(3)-Vend2(2)]% case alpha=1.3

abs(v2(1)-v2(2))


v3=[2*Vend3(2)-Vend3(1), 2*Vend3(3)-Vend3(2)]% case alpha=1.5

abs(v3(1)-v3(2))


v4=[2*Vend4(2)-Vend4(1), 2*Vend4(3)-Vend4(2)]% case alpha=1.7

abs(v4(1)-v4(2))
