simpsons rule for integration in matlab

main.m
function I = simpsons_rule(f, a, b, N)
% SIMPSONS_RULE Calculates definite integral of f(x) from a to b 
% using Simpson's 1/3 rule with N intervals.
% 
% Inputs:
%   - f: function handle for the function to be integrated
%   - a: lower limit of integration
%   - b: upper limit of integration
%   - N: number of intervals (must be even)
%
% Output:
%   - I: definite integral value

% Check if N is even:
if mod(N, 2) ~= 0
    error('N must be even for Simpson''s 1/3 rule.')
end

% Define interval width:
h = (b - a) / N;

% Calculate sum of odd terms and even terms:
odd_sum = 0;
even_sum = 0;

for i = 1:N-1
    x_i = a + i*h;
    if mod(i, 2) == 0
        even_sum = even_sum + f(x_i);
    else
        odd_sum = odd_sum + f(x_i);
    end
end

% Calculate integral value:
I = (h/3) * (f(a) + 4*odd_sum + 2*even_sum + f(b));

end
834 chars
39 lines

main.m
f = @(x) x^2;
a = 0;
b = 2;
N = 4;

I = simpsons_rule(f, a, b, N);
67 chars
7 lines