composite simpsons method in matlab

main.m
function I = composite_simpson(f,a,b,n)
% This function performs composite Simpson's method for numerical integration.
% Inputs:
%   f - the integrand as a function handle
%   a - the lower bound of integration
%   b - the upper bound of integration
%   n - the number of subintervals (must be even)
% Output:
%   I - the approximated value of the integral
h = (b-a)/n; % width of each subinterval
x = a:h:b; % vector of x values
y = f(x); % vector of corresponding y values
% Simpson's rule formula:
I = h/3*(y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2)) + y(end));
% Composite Simpson's rule formula:
for i = 1:n/2-1
    a_i = a + 2*i*h; % lower bound of ith subinterval
    b_i = a + (2*i+2)*h; % upper bound of ith subinterval
    x_i = a_i:h:b_i; % vector of x values for ith subinterval
    y_i = f(x_i); % vector of corresponding y values for ith subinterval
    I = I + h/3*(y_i(1) + 4*sum(y_i(2:2:end-1)) + 2*sum(y_i(3:2:end-2)) + y_i(end)); % add contribution of ith subinterval
end
996 chars
23 lines

main.m
f = @(x) 1./x.^2; % example integrand
a = 1; % example lower bound of integration
b = 2; % example upper bound of integration
n = 10; % example number of subintervals
I = composite_simpson(f,a,b,n); % call composite Simpson's method function
disp(I); % display approximated value of integral
292 chars
7 lines