function res = simpson(f,a,b,tol)
%  res = simpson(fun,a,b,tol) evaulate an integral by 
%    adaptive Simpson's method, using the algorithm on page 241.
%    A plot illustrates the adaptive process.
%    The partitioning is given by the red points. 

level_max = 10;
hold off
res = simpson_plot(f,a,b,tol,0,10)
hold on

%---------------------------------------------------------------
function simpson_result = simpson_plot(f,a,b,tol,level,level_max) 
if level==0
    fplot(f,[a,b]);
    grid on
    hold on
    pause(1)
end
level = level+1
h = b-a;
c = (a+b)/2;
one_simpson = h*(feval(f,a)+4*feval(f,c)+feval(f,b))/6;
d = (a+c)/2;
e = (c+b)/2;
plot(a,feval(f,a),'*r',b,feval(f,b),'*r');
plot(c,feval(f,c),'ob',d,feval(f,d),'og',e,feval(f,e),'og');
two_simpson = h*(feval(f,a)+4*feval(f,d)+2*feval(f,c)+4*feval(f,e) ...
		 + feval(f,b))/12;
errtemp = abs(two_simpson-one_simpson)/15;
title(sprintf('tol=%5.2e, error=%5.2e', tol, errtemp));
pause(1)
if level>=level_max
    simpson_result = two_simpson;
    warning('Maximum level reached')
else
    err = (two_simpson-one_simpson)/15;
    if abs(err)<tol 
	simpson_result = two_simpson+(two_simpson-one_simpson)/15;
    else
	left_simpson = simpson_plot(f,a,c,tol/2,level,level_max);
	right_simpson = simpson_plot(f,c,b,tol/2,level,level_max);
	simpson_result = left_simpson+right_simpson;
    end
end

