function for simpns in matlab

function [I] = simpson(f,a,b,n)
% Simpson's 1/3 rule for numerical integration

% Input:
% f - the integrand
% a - the lower limit of integration
% b - the upper limit of integration
% n - the number of intervals to use (must be even)

% Output:
% I - the approximate value of the integral

% Check that n is even:
if (mod(n,2) ~= 0)
    error('The number of intervals (n) must be even');
end

% Calculate the interval width and x values:
h = (b - a) / n;
x = a:h:b;

% Calculate the function values:
y = f(x);

% Apply 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) );
end
621 chars
28 lines

% Example usage:
f = @(x) x.^2 .* exp(-x); % the integrand
a = 0; % lower limit of integration
b = 2; % upper limit of integration
n = 8; % number of intervals (must be even)
I = simpson(f,a,b,n) % approximate value of the integral
232 chars
7 lines