simpsons rule using recursion in matlab

main.m
function I = SimpsonsRule(f,a,b,n)
% This function uses Simpson's Rule to numerically approximate the integral of
% f(x) over the interval [a,b]. It recursively divides the interval [a,b] into
% n subintervals and applies Simpson's Rule to each subinterval.
%
% Input:
%   f: the function to be integrated (must be a vectorized function)
%   a,b: the limits of integration
%   n: the number of subintervals to use (must be an even integer)
%
% Output:
%   I: the numerical approximation of the integral of f(x) over [a,b]

h = (b-a)/n; % width of subintervals
x = a:h:b; % vector of evaluation points

if length(x) ~= n+1
    error('n must be an even integer.')
end

fx = f(x); % evaluate f(x) at each point in x

% evaluate integral using Simpson's Rule
I = (h/3) * (fx(1) + 4*sum(fx(2:2:end-1)) + 2*sum(fx(3:2:end-2)) + fx(end));

% check for convergence
if n == 2
    return
else
    % divide interval into two subintervals
    n = n/2;
    left_I = SimpsonsRule(f,a,a+h,n);
    right_I = SimpsonsRule(f,a+h,b,n);
    I = left_I + right_I;
end

end
1052 chars
38 lines