File: ef230-2023-08/www/examples/simpson.m Download
%! Approximate integral using composite simpson rule
%! http://en.wikipedia.org/wiki/Simpson's_rule
function r=simpson(f,a,b,n)
% f - handle of function to be integrated
% a - lower bound
% b - upper bound
% n - number of intervals, defaults to 100, must be even
% r - result of integration

if nargin<4
    n=100;
end
if nargin<3
    error('insufficient inputs');
end
if mod(n,2)==1
    error('n must be even');
end
if a>b
    error('b must be greater than a');
end
h=(b-a)/n; % interval size
x=linspace(a,b,n+1); % n+1 points
y=f(x); % calc y values
j=1:(n/2-1); % j according to given formula
k=1:(n/2); % k according to given formula
% Composite Simpson Rule
% the extra +1 accounts for Matlab start index of 1, not 0
r=(h/3).*(f(a) + f(b) + 2*sum(y(2*j+1)) + 4*sum(y(2*k-1+1)));
return