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

% We will integrate over the interval from 0 to 1
n    = 101; % 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;

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

% One-dimensional construction is completed!
% Now we need to use meshgrid to go to 2D
[X,Y]   = meshgrid(xvec);
[W1,W2] = meshgrid(w);
Wmat    = W1.*W2;

% Construt the function to be integrated
fun = sqrt(X).*Y;
%fun = X.^2.+Y.^2;
%fun = ones(size(X));
fun = fun(:);

% Integrate by inner product
int = (Wmat(:).')*fun

% Take a look at the points and the function
figure(1)
clf
subplot(1,2,1)
plot(X(:),Y(:),'r.')
axis square
subplot(1,2,2)
mesh(X,Y,reshape(fun,size(X)))
axis square











