%MATLAB version of the script for example 7.9 on pages 343f of conte/deboor
% example 7.9  Romberg integration.
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
F = @(x) exp(-x.*x);
disp('Example 7.9: Romberg T-table')

A = 0; B = 1;  MSTART = 2; NROW = 6;
T = ROMBRG( F, A, B, MSTART, NROW );

function T = ROMBRG( F, A, B, MSTART, NROW )

M = MSTART; T = zeros(NROW,NROW);
T(1,1) = TRAP(F,A,B,M);
disp('  Romberg T-table')
form = [ ' %15.7e ', ' %15.7e ', ' %15.7e ', ' %15.7e ', ' %15.7e ', ...
         ' %15.7e ', ' %15.7e '];
fprintf([form(1:8), '\n'], T(1,1))
for k=2:NROW
   T(k,1) = (T(k-1,1)+MIDPOINT(F,A,B,M))/2; M = 2*M;
     % save differences for later calc. of ratios
   for j=1:k-1
      T(k-1,j) = T(k,j) - T(k-1,j);
      T(k,j+1) = T(k,j) + T(k-1,j)/(4^j-1);
   end
   fprintf([form(1:8*k),'\n'], T(k,1:k))
end
    % calculate ratios
disp(' Table of Ratios')
for k=1:NROW-2
   index = T(k+1,1:k) == 0;
   T(k,~index) = T(k,~index)./T(k+1,~index); 
   T(k,index) = 0;
   fprintf([form(1:8*k),'\n'], T(k,1:k))
end
return, end     


function value = MIDPOINT( F, A, B, N )
H = (B-A)/N; value = H*sum(F(A+H*((0:N-1)+1/2)));
return, end
