function progex3
    fprintf('\na)\n\n');
    maketable('4 / (1 + x^2)', 0, 1, 10, 1e-12, pi);

    fprintf('\nb)\n\n');
    maketable('4*x^2*exp(-x^2)', 0, 7, 10, 1e-12, sqrt(pi));

    fprintf('\nc)\n\n');
    maketable('4*x^2*exp(-x^2)', 0, 500, 12, 0, sqrt(pi));

    fprintf('\nd)\n\n');
    maketable('sqrt(x*(4-x))', 0, 2, 10, 1e-12, pi);

    fprintf('\ne)\n\n');
    maketable('sin(64*x)^2', 0, 2*pi, 12, 0, pi);

    fprintf('\nf)\n\n');
    maketable('1/sqrt(1 - sin(x)^2/2)', 0, pi, 10, 1e-12, 3.7081493546027);

    fprintf('\ng)\n\n');
    maketable('1/sqrt(1 - sin(x)^2/2)', 0, 1, 10, 1e-12, 1.0832167728452);
end

function maketable(F, a, b, n, Tol, I)
    f = inline(vectorize(F));
    
    T = richardson(f, a, b, n, Tol);
    
    n = size(T,1);
    
    T(:,5) = abs(I - T(:,1));
    T(2:n,6) = abs(I - T(2:n,2));
    
    S = zeros(n, 4) + nan;
    S(4:n,1) = abs(T(4:n,4) - T(4:n,3));
    S(4:n,2) = (T(3:n-1,2) - T(2:n-2,2)) ./ (T(4:n,2) - T(3:n-1,2));
    S(5:n,3) = (T(4:n-1,3) - T(3:n-2,3)) ./ (T(5:n,3) - T(4:n-1,3));
    S(6:n,4) = (T(5:n-1,4) - T(4:n-2,4)) ./ (T(6:n,4) - T(5:n-1,4));
    
    fprintf('          T0          T1          T2          T3    |I - T0|    |I - T1|\n');
    fprintf('%12.3e%12.3e%12.3e%12.3e%12.3e%12.3e\n', T');
    fprintf('\n   |T3 - T2|     Ratio 1     Ratio 2     Ratio 3\n');
    fprintf('%12.3e%12.3e%12.3e%12.3e\n', S');
end

function T = richardson(f, a, b, n, Tol)

    % compute trapezoid rule with 1 subdivision

    T = zeros(n, 4) + nan;
    T(1, 1) = (f(a) + f(b))*(b - a)/2;
    
    for i = 2:n
        % efficiently compute next trapezoid entry
        
        h = (b - a)/2^(i-1);
        x = a + h*(1:2:2^(i-1)-1);
        T(i, 1) = T(i-1, 1)/2 + h*sum(f(x));
        
        % compute all the first 3 extrapolants
        
        for j = 2:min(i, 4)
           T(i, j) = (4^(j-1) * T(i, j-1) - T(i-1, j-1)) / (4^(j-1) - 1);
        end

        % check Tol
        
        if i > 4 && Tol ~= 0 && abs(T(i,4) - T(i-1,4)) < Tol
            T = T(1:i,:);
            break;
        end
    end
    
end