function simpson_result = simpson(f,a,b,tol,level,level_max)
%
% simpson_result = simpson(fun,a,b,tol,0,level_max)
%
% Regner ut integralet av fun fra a til b 
% med adaptiv Simpson's metode. 
%
% level_max=10 er en passende verdi. 
%
level = level+1;
h = b-a;
c = (a+b)/2;
S1 = h*(feval(f,a)+4*feval(f,c)+feval(f,b))/6;
d = (a+c)/2;
e = (c+b)/2;
S2 = h*(feval(f,a)+4*feval(f,d)+2*feval(f,c)+4*feval(f,e) ...
		 + feval(f,b))/12;
if level>=level_max
    simpson_result = S2;
    warning('Maximum level reached');
else
    err = (S2-S1)/15;
    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

