% Demonstration of Simpson's rule
%
% Samuli Siltanen November 2012

% We will integrate over the interval from 0 to 1
n    = 21; % Must be odd!
xvec = [0:(n-1)]*(1/(n-1));

% Construct Bart Simpson's hair
tmp1 = 4*ones(1,(n-1)/2);
tmp2 = 2*ones(1,(n-1)/2);
w    = [tmp1;tmp2];
w    = [1;w(:)];
w(end) = 1;

% Construt the function to be integrated
%fun = xvec.^2;
fun = sqrt(xvec);
fun = fun(:);

% Construct quadrature weights
w   = w/sum(w);
w   = w(:).';

% Integrate by inner product
int = w*fun



