% Demonstration of Simpson's rule
% We will integrate over the interval from 0 to 1
%
% Samuli Siltanen November 2012

% Use this vector to store the values of the integral at different levels
% of approximation. Warning: this method should not be used for big
% matrices!
intvec = [];

% Vector of n values. The numbers n must be odd!
nvec = 11:10:521;

% Loop over approximations to study convergence
for n = nvec
    
    % Construct the quadrature points
    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
    intvec = [intvec;w*fun];
    
end

% Display the results
format long
disp([nvec(:),intvec])
