Adaptive Simpson's rule for a bivariate function

adapsimpson2.m

clear;

T = 7;
f = @(x, y) exp(-T * (x.^2 - y.^2)) - exp(-T * (x.^2 + y.^2));
xa = 0;
xb = 2;
ya = 0;
yb = 3;

S = adapsimpson2(f, xa, xb, ya, yb, 1.e18)

function S = simpson(f, x0, x2)
    h = (x2 - x0) / 2;
    x1 = x0 + h;
    S = (h / 3) * (f(x0) + 4 * f(x1) + f(x2));
end

function S = adapsimpson(f, x0, x2, tol)
    h = (x2 - x0) / 2;
    x1 = x0 + h;
    S1 = simpson(f, x0, x2);
    S2 = simpson(f, x0, x1) + simpson(f, x1, x2);
    E2 = (S2 - S1) / 15;
    if abs(E2) < tol
        S = S2 + E2;
    else
        Sl = adapsimpson(f, x0, x1, tol);
        Sr = adapsimpson(f, x1, x2, tol);
        S = Sl + Sr;
    end
end

function S = adapsimpson2(f, xa, xb, ya, yb, tol)
    % g is the integral of f from ya to yb with x fixed, solved numerically
    % with adaptive Simpson's.
    % The value of x is determined by adapsimpson, which evaluates g(x0),
    % g(x1) and g(x2) for a given interval [x0, x2], initially [xa, xb].
    g = @(x) adapsimpson(@(y) f(x, y), ya, yb, tol);
    S = adapsimpson(g, xa, xb, tol);
end