adaptive simpson matrix multiplication in matlab

adaptive_simpson.m

function [S, count] = adaptive_simpson(f, a, b, epsilon, level, level_max)
% Evaluates the integral of f using Adaptive Simpson
% over [a, b] such that 1 / 15 * abs(S^(2) - S^(1)) < e.
count = 1;
level = level + 1;
h = b - a;
c = (a + b) / 2;
f_a = f(a);
f_c = f(c);
f_b = f(b);
one_simpson = h * (f_a + 4 * f_c + f_b) / 6;
d = (a + c) / 2;
e = (c + b) / 2;
f_d = f(d);
f_e = f(e);
two_simpson = h * (f_a + 4 * f_d + 2 * f_c + 4 * f_e + f_b) / 12;
if level >= level_max
    S = two_simpson;
    disp('maximum level reached');
else
    if abs(two_simpson - one_simpson) < (15 * epsilon)
        S = two_simpson + (two_simpson - one_simpson) / 15;
    else
        [left_simpson, left_count] = adaptive_simpson(f, a, c, epsilon / 2, level, level_max);
        [right_simpson, right_count] = adaptive_simpson(f, c, b, epsilon / 2, level, level_max);
        S = left_simpson + right_simpson;
        count = left_count + right_count;
    end
end
end

adaptive_simpson_2d.m

function S = adaptive_simpson_2d(f, a, b, c, d, epsilon, level, level_max)
% Evaluates the double integral of f(x, y) using Adaptive Simpson
% over [a, b] such that 1 / 15 * abs(S^(2) - S^(1)) < e.
syms y;
level = level + 1;
h = b - a;
i = (a + b) / 2;
f_a = @(y) f(a, y);
f_i = @(y) f(i, y);
f_b = @(y) f(b, y);
F_a = adaptive_simpson(f_a, c, d, epsilon, level, level_max);
F_i = adaptive_simpson(f_i, c, d, epsilon, level, level_max);
F_b = adaptive_simpson(f_b, c, d, epsilon, level, level_max);
one_simpson = h * (F_a + 4 * F_i + F_b) / 6;
j = (a + i) / 2;
k = (i + b) / 2;
f_j = @(y) f(j, y);
f_k = @(y) f(k, y);
F_j = adaptive_simpson(f_j, c, d, epsilon, level, level_max);
F_k = adaptive_simpson(f_k, c, d, epsilon, level, level_max);
two_simpson = h * (F_a + 4 * F_j + 2 * F_i + 4 * F_k + F_b) / 12;
if level >= level_max
    S = two_simpson;
    disp('maximum level reached');
else
    if abs(two_simpson - one_simpson) < (15 * epsilon)
        S = two_simpson + (two_simpson - one_simpson) / 15;
    else
        left_simpson = adaptive_simpson_2d(f, a, i, c, d, epsilon / 2, level, level_max);
        right_simpson = adaptive_simpson_2d(f, i, b, c, d, epsilon / 2, level, level_max);
        S = left_simpson + right_simpson;
    end
end
end