% 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 = 10:10:100;

% Loop over approximations to study convergence
for n = nvec
    
    % Construct the quadrature points and weights
    [xvec,w] = gaussint(n,0,1);
    w        = w(:).';
    xvec     = xvec(:);
        
    % Construt the function to be integrated
    %fun = xvec.^2;
    fun = sqrt(xvec);
    fun = fun(:);
    
    % Integrate by inner product
    intvec = [intvec;w*fun];
    
end

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