simpson's 1/3 rule in matlab in matlab

Here's an implementation of Simpson's 1/3 rule in MATLAB:

main.m
function I = simpsons13(fun, a, b, n)
% Calculates the definite integral of a function using Simpson's 1/3 rule
% fun: function handle to the integrand
% a: lower limit of integration
% b: upper limit of integration
% n: number of subintervals (must be even)
  
  if mod(n,2) ~= 0
    error('n must be an even integer');
  end
  
  x = linspace(a,b,n+1);
  y = fun(x);
  h = (b-a)/n;
  
  % Applying Simpson's 1/3 rule
  I = h/3 * (y(1) + 4*sum(y(2:2:end)) + 2*sum(y(3:2:end-2)) + y(n+1));
end
494 chars
19 lines

This function takes four inputs: a function handle to the integrand fun, the lower and upper limits of integration a and b, and the number of subintervals n. It first checks if n is even, and throws an error otherwise.

Then, it creates the array of points x that split the interval [a,b] into n subintervals (with n+1 points). It evaluates the integrand at these points to obtain the array y.

Finally, it applies Simpson's 1/3 rule to compute the approximation for the integral. The formula is applied directly using array operations, with special treatment for the first and last points.

Here's an example of how to use this function:

main.m
% Define the integrand
f = @(x) exp(-x.^2);

% Compute the integral from 0 to 1 using 10 subintervals
I = simpsons13(f, 0, 1, 10);
131 chars
6 lines