function [I,x,y] = adSimpson(f,a,b,TOL,level,level_max)
% Approximates an integral using adaptive Simpson's quadrature
%
% INPUT: f (function to be integrated)
%        a (left endpoint)
%        b (right endpoint)
%        TOL (error tolerance)
%        level (recursion level of current interval, starts at 0)
%        level_max (maximum recursion level depth)
%
% OUTPUT: I (Integral approximation)
%         x (Vector of partition points.)
%         y (Vector of function values at partition points)
%
level = level + 1;
h = b-a;
c = (a+b)/2;
S1 = (h/6)*(f(a)+4*f(c)+f(b));

d = (a+c)/2;
e = (c+b)/2;

S2 = (h/12)*(f(a)+4*f(d)+2*f(c)+4*f(e)+f(b));

if (1/15)*abs(S2-S1) < TOL
    I = S2 + (S2-S1)/15;
    x = [a,d,c,e,b];
    y = f(x);
    
elseif level >= level_max
    I = S2;
    fprintf('Max level reached\n\n');
    x = [a,d,c,e,b];
    y = f(x);
    
else
    [L,xL,yL] = adSimpson(f,a,c,TOL/2,level,level_max);
    [R,xR,yR] = adSimpson(f,c,b,TOL/2,level,level_max);
    I = L+R;
    x = [xL,xR(2:end)];
    y = [yL,yR(2:end)];
end


end

