Adaptive numerical integration

Example for integral

We want to approximate the integral
with an approximation Q such that Q-I<=1e-2. We want to use the composite trapezoid rule with a subdivision which is adapted to the behavior of the function f: large intervals where |f''| is small, small intervals where |f''| is large.
clearvars;
global Nev % use global variable Nev to count function evaluations
f = @(x) x./(1+x.^4);
a = 0; b = 4;
exact = (atan(b^2)-atan(a^2))/2;
fplot(f,[a,b])
xlabel('x'); title('function f(x)')

Try adaptint with tol=.1, .01, .01

Note that the actual error is about 0.3 times smaller than the tolerance.
Try to change the line Q=QT to Q=QS in the code below. What happens with the errors?
for tol = [.1 .01 .001]
figure; fplot(f,[a,b]); hold on
Nev = 2; % 2 initial evaluations at a,b
Q = adaptint(f,a,b,tol)
Nev % number of function evaluations
error = Q - exact
title(sprintf('tol = %g: adaptint uses %g subintervals',tol,Nev-1))
hold off;
end
Q = 0.79274
Nev = 13
error = 0.038549
Q = 0.75704
Nev = 41
error = 0.0028541
Q = 0.75439
Nev = 113
error = 0.00020184

Recursive function adaptint(f,a,b,Tol)

We use the composite trapezoid rule. On each subinterval we estimate the error |QT-I| by |QT-QS| where QS is the Simpson rule value. We subdivide the intervals until the sum of estimated errors is <= Tol.
In order to demonstrate how this works we
function Q = adaptint(f,a,b,Tol)
global Nev
fa = f(a); fb = f(b);
QT = (b-a)/2*(fa+fb); % Trapezoid rule
c = (a+b)/2; fc = f(c); % evaluate f in midpoint
Nev = Nev + 1; % count function evaluations
plot([c c],[0,2]*1e-2,'r'); % plot new point
QS = (b-a)/6*(fa+4*fc+fb); % Simpson rule
% for small intervals QS is much closer to I than QT, hence QT-I ~ QT-QS
if abs(QT-QS)<=Tol % if estimated error is <= Tol
Q = QT; % accept trapezoid rule value
else
Q = adaptint(f,a,c,Tol/2) + adaptint(f,c,b,Tol/2); % use algorithm for [a,c] and [c,b]
end
end