function [q, x] = simpson_adat (fun, a, b, tol)
%
% function [q, x] = simpson_adat (fun, a, b, tol)
%
% Formula di Simpson adattativa in [a,b].
q = simpson (fun, 2, a, b);
[q1, x1] = simpson (fun, 2, a, (a + b) / 2);
[q2, x2] = simpson (fun, 2, (a + b) / 2, b);
if (abs (q1 + q2 - q) > 10 * tol)
   [q1, x1] = simpson_adat (fun, a, (a + b) / 2, tol / 2);
   [q2, x2] = simpson_adat (fun, (a + b) / 2, b, tol / 2);
end
q = q1 + q2;
x = [x1; x2];
%!demo
%! fun = @(x) 100 ./ x .^ 2 .* sin (10 ./ x);
%! a = 1;
%! b = 3;
%! [q, x] = simpson_adat (fun, a, b, 1e-4);
%! plot(x, fun (x), x, zeros (size (x)), 'o')
%! E = abs (q - quad (fun, a, b))
