create a function using repeated simpsons rule in matlab

main.m
function I = repeated_simpsons_rule(f, a, b, n)
%Repeated Simpsons rule for numerical integration
%Inputs:
%   f - function to be integrated
%   a - integration interval start
%   b - integration interval end
%   n - number of intervals to be used
%Output:
%   I - numerical estimate of the integral of f over [a,b]

if mod(n,2) == 1
    error('n must be an even positive integer')
end

h = (b-a)/n;
x = a:h:b;
y = f(x);

%Repeating Simpson's rule for each interval
I = y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2))...
    + 4*sum(y(4:2:end-4)) + y(end);

I = I*h/3;
end
574 chars
25 lines