function ex17()
fB = @(x) cos(10*exp(-2*(x-1).^2));
a = 0.4;
b = 4.0;
tol = 0.01;

fprintf("   a      b       I1      I2\n")
fprintf("===========================================\n")
[I points] = simpson_adapt(fB, a, b, tol);
Igenau = quad(fB, a, b, 1e-8);

fprintf("\nIntegral: %12.8f  genau: %12.8f\n", I, Igenau);
fprintf("%d Punkte: \n", length(points));
fprintf("   %5.3f %5.3f %5.3f %5.3f\n", points);
fprintf("\n");

x = a:0.001:b;
y = fB(x);
plot(x, y, "k-", x, 0*x, "k-", points, 0*points, "r*")
axis([a b -1.1 1.1]);

F = getframe(gcf);
imwrite(F.cdata, "bild60.png");

%--------------------------------------------------
function [I points] = simpson_adapt(f, a, b, tol)
% berechnet das Integral mit adaptivem Simpsonverfahren
% Eingabe
%      f   Integrand
%      a   untere Intervallgrenze
%      b   obere Intervallgrenze
%    tol   Toleranz
% Ausgabe
%      I   Wert des Integrals
%  points  Liste der Zwischenpunkte

I1 = (b-a)/6 * (f(a) + 4*f((a+b)/2) + f(b));
I2 = (b-a)/12 * (f(a) + 4*f((3*a+b)/4) + 2*f((a+b)/2) ...
                  + 4*f((a+3*b)/4) + f(b));
fprintf("%6.3f %6.3f  %7.4f %7.4f", a, b, I1, I2)

if  abs(I1 - I2) < tol
  fprintf("   ok\n");
  I = (16*I2 - I1)/15;
  points = [a b];
else
  fprintf("   halbieren\n");
  c = (a+b)/2;
  [Il pointsl] = simpson_adapt(f, a, c, tol);
  [Ir pointsr] = simpson_adapt(f, c, b, tol);
  I = Il + Ir;
  points = [pointsl pointsr(2:end)];
end
