function simpson_result = simpson(f,a,b,tol,level,level_max)
%
% simpson_result = simpson(f,a,b,tol,0,level_max)
%
% Compute the integral of a function f from a to be within a tolerance tol,
% using an adaptiv Simpson method. 
%
% level_max=10 is a suitable value. 
%
level = level+1;
h = b-a;
c = (a+b)/2;
S1 = h*(f(a)+4*f(c)+f(b))/6;
d = (a+c)/2;
e = (c+b)/2;
S2 = h*(f(a)+4*f(d)+2*feval(f,c)+4*f(e) + f(b))/12;
if level>=level_max
    simpson_result = S2;
    warning('Maximum level reached');
else
    err = (S2-S1)/15;
    a,b, S1, S2, err, pause
    if abs(err)<tol 
        simpson_result = S2+err
    else
        left_simpson = simpson(f,a,c,tol/2,level,level_max);
        right_simpson = simpson(f,c,b,tol/2,level,level_max);
        simpson_result = left_simpson+right_simpson;
    end
end

