clear
format short e
format compact
%%
% We estimate \int_0^2 x^2 e^{-2x} dx using extrapolation. 
f = @(x) x.^2.*exp(-2*x);
a = 0;  
b = 2;  
I = integral(f, a, b, 'abstol', 1e-14, 'reltol', 1e-14);

%%
% We start with the trapezoid formula on n nodes.
n = 20;  
[T(1), ~, ~] = trapezoid(f, a, b, n);
err_2nd = abs(I - T(1))

%%
% Now we double to n=2N
n = 2*n;  
[T(2), ~, ~] = trapezoid(f, a, b, n);
err_2nd = abs(I - T(2))

%%
% As expected for a second-order estimate, the error went down by a factor
% of about 4. We can repeat the same code to double n again.
n = 2*n;  
[T(3), ~, ~] = trapezoid(f, a, b, n);
err_2nd = abs(I - T(3))

%%
% Let us now do the first level of extrapolation to get results from
% Simpson's formula. We combine the elements |T(i)| and |T(i+1)| the same
% way for i=1 and i=2.
S = (4*T(2:3) - T(1:2)) / 3;
err_4th = abs(I - S)

%%
% With the two Simpson values S(N) and S(2N) in hand, we can do one
% more level of extrapolation to get a 6th-order accurate result.
R = (16*S(2) - S(1)) / 15;
err_6th = abs(I - R)

%%
% If we consider the computational time to be dominated by evaluations of
% f, then we have obtained a result with twice as many accurate digits as
% the best trapezoid result, at virtually no extra cost.