function I = Romberg( f, a, b, n )
%Romberg:  Approximate integral via Romberg integration
% 
% >> I = Romberg( f, a, b, n )
% 
% in:  f    = function to be integrated (name or handle)
%      a, b = limits of integration
%      n    = number of extrapolations to do
% 
% out: I    = final, fully-extrapolated approximation to integral
% 
% output to terminal: T, E, and R tables (extrapolations, error estimates,
%                                         and ratios)
% 
% examples:  >> I = Romberg( 'sin', 0, pi/2, 5 )
%            >> I = Romberg( @sin , 0, pi/2, 5 )
%            >> I = Romberg( inline('sin(x)'), 0, pi/2, 5 )

% ~/matlab/courses/NumerComput/Romberg.m
% 
% C.G. (3-02)

% build vector of 'n' successive composite-Trapezoid-rule approximations

ncells = 1 ;  h = b - a ;              % initialize # of cells and cell width

x = [ a, b ] ;                         % initial grid = 1 cell, whole interval

T(1) = h * sum( feval( f, x ) ) / 2 ;  % basic Trap rule over whole interval

for i = 1 : n

  ncells = 2 * ncells ;  h = h / 2 ;            % refined # cells and width

  xc = ( x(1:end-1) + x(2:end) ) / 2 ;          % cell centers of current grid

  T(i+1) = T(i) / 2 + h * sum( feval(f,xc) ) ;  % next comp-Trap value

  npts = ncells + 1 ;

  x(1:2:npts  ) = x ;                           % refined grid
  x(2:2:npts-1) = xc ;

end

% extrapolate

[ Textrap, T, E, R ] = extrap( T ) ;

I = Textrap ;                          % return fully extrapolated value

% print extrapolation tables to terminal window

fprintf( 1, '\n' )                              % T-table
fprintf( 1, 'T-table \n' )
fprintf( 1, '------- \n' )

for i = 1 : size(T,1)
  fprintf( 1, ' %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f ', T(i,1:i) )
  fprintf( 1, '\n' )
end

fprintf( 1, '\n' )
fprintf( 1, 'E-table \n' )                      % E-table
fprintf( 1, '------- \n' )

for i = 1 : size(E,1)
  fprintf( 1, ' %9.2e %12.2e %12.2e %12.2e %12.2e ', E(i,1:i) )
  fprintf( 1, '\n' )
end

fprintf( 1, '\n' )
fprintf( 1, 'R-table \n' )                      % R-table
fprintf( 1, '------- \n' )

for i = 1 : size(R,1)
  fprintf( 1, ' %5.2f %12.2f %12.2f %12.2f ', R(i,1:i) )
  fprintf( 1, '\n' )
end

%-----------------------------------------------------------------------------

function [ Textrap, T, E, R ] = extrap( Tvec )
%extrap:  extrapolate a sequence to zero mesh width ... assuming bisection
%         refinement and an asymptotic expansion including only EVEN powers:
% 
%           T(h) = T(0) + c1*h^2 + c2*h^4 + c3*h^6 + ...
% 
% usage:  [ Textrap, T, E, R ] = extrap( Tvec )
% 
% in:   Tvec    = vector of values to extrapolate
% 
% out:  Textrap = fully extrapolated value
%       T, E, R = the extrapolation, error, and ratio tables

I = length( Tvec ) ;

T = zeros(I,I) ;  E = zeros(I-1,I-1) ;  R = zeros(I-2,I-2) ; % allocate arrays

T(:,1) = Tvec(:) ;           % initialize 1st col of T-array with input values

% build the T and E tables column by column

for j = 1 : I-1
  
  i = j : I-1 ;
  
  E(i,j) = ( T(i+1,j) - T(i,j) ) / ( 4^j - 1 ) ;    % est of error in T(i+1,j)
  
  T(i+1,j+1) = T(i+1,j) + E(i,j) ;                  % corrected/improved value
  
end

Textrap = T(I,I) ;

% build the R table column by column ... gives the ratios of successive
% error estimates ... if column is converging like O(h^p), then this 
% ratio should be 2^p

for j = 1 : I-2
  
  i = j : I-2 ;
  
  R(i,j) = E(i,j) ./ E(i+1,j) ;
  
end