quadrature integration using recursion in matlab

main.m
function integral = quad_integration(f, a, b, tol)
% QUAD_INTEGRATION Computes the quadrature integration of a function
% using recursive sub-division
% Usage:
%   integral = quad_integration(f, a, b, tol)
%
%   f - Function to integrate
%   a - Lower bound of the integration interval
%   b - Upper bound of the integration interval
%   tol - Tolerance for accuracy of the result

% Calculate the midpoint of the interval
mid = (a+b)/2;

% Calculate the function values at each endpoint and the midpoint
fa = f(a);
fb = f(b);
fm = f(mid);

% Calculate the trapezoidal integration of the whole interval
trap = (b-a)*(fa + fb)/2;

% Calculate the Simpson's integration of the whole interval
simp = (b-a)*(fa + 4*fm + fb)/6;

% Calculate the absolute difference between Simpson's and trapezoidal
diff = abs(simp - trap);

% If the difference is within the given tolerance, use Simpson's rule
if diff <= tol
    integral = simp;

% Otherwise, recursively sub-divide the interval and repeat
else
    integral = quad_integration(f, a, mid, tol/2) ...
             + quad_integration(f, mid, b, tol/2);
end
1101 chars
38 lines