% cs3734s04HW09.m
% composite quadrature rules
% K. Ming Leung, 05/03/04

clear all; format long e;
%fcn = 'cos'; a = 0; b = pi/2; Iexact = 1;
fcn = inline('4./(1+x.*x)'); %'Lorenzian'; 
a = 0; b = 1; Iexact = pi;
ba = b-a;
n = 1;
fprintf('\n');
fprintf('         h        Midpoint Error  Trapezoid Error  Simpson Error\n');
fprintf('================================================================\n');

for loop = 1:19
    n = 2*n; h = ba/n;
    Error(loop,1) = h;
    Imidpoint = midpoint(fcn,a,b,n);
    Error(loop,2) = Imidpoint-Iexact;

    Itrapezoid = trapezoid(fcn,a,b,n);
    Error(loop,3)  = Itrapezoid-Iexact;

    ISimpson = Simpson(fcn,a,b,n);
    Error(loop,4)  = ISimpson-Iexact;
    fprintf('%15.4e %15.4e %15.4e %15.4e \n',Error(loop,:));
end

% Romberg quadrature
IRomberg = Romberg(fcn,a,b,7),  % 7 iterations for the best result
errorRomb = IRomberg(end,end)-Iexact,
