% demonstration of calling romberg.m 
% based on integrating sin(x) from 0 to pi, 
% which has the exact answer 2.

 a     = 0; 
 b     = pi;
 max_k = 6;

 R = romberg('sin',a,b,max_k);

% The ratio of the difference between successive entries in column k 
% should be approximately 4^k.  This provides a simple test to see if
% our extrapolation is working, without requiring knowledge of the exact
% answer.  See Johnson and Riess, Numerical Analysis, 2nd ed., p. 323.

 Ratios = zeros(max_k-1,max_k-1);
 for k=1:max_k-1
    Ratios(k:end,k) = (R(1+k:end-1,k)-R(k:end-2,k))./(R(2+k:end,k)-R(1+k:end-1,k));
 end

 format long
 R
 Ratios

 break
 fprintf('%13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f\n', ...
         R(1,1),R(1,2),R(1,3),R(1,4),R(1,5),R(1,6),R(1,7));
 fprintf('%13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f\n', ...
         R(2,1),R(2,2),R(2,3),R(2,4),R(2,5),R(2,6),R(2,7));
 fprintf('%13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f\n', ...
         R(3,1),R(3,2),R(3,3),R(3,4),R(3,5),R(3,6),R(3,7));
 fprintf('%13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f\n', ...
         R(4,1),R(4,2),R(4,3),R(4,4),R(4,5),R(4,6),R(4,7));
 fprintf('%13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f\n', ...
         R(5,1),R(5,2),R(5,3),R(5,4),R(5,5),R(5,6),R(5,7));
 fprintf('%13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f\n', ...
         R(6,1),R(6,2),R(6,3),R(6,4),R(6,5),R(6,6),R(6,7));
 fprintf('%13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f  %13.12f\n', ...
         R(7,1),R(7,2),R(7,3),R(7,4),R(7,5),R(7,6),R(7,7));