function useromberg %%%%% test 1
   format long e
   format compact
   n=6
   f = @exp
   a = 0, b=1
   tableau=romberg(f,a,b,n)
   exact = exp(1)-exp(0)
   exactm = tril(exact*ones(n,n));
   relerr = (tableau-exactm)/exact
end

function useromberg2   %%%%% test2
   format short e
   format compact
   n=8
   f = @(x) 1/(1+12*x^2)    % Runge's function
   a = 0, b = 1
   exact = atan(2*sqrt(3))*sqrt(3)/6
   tableau=romberg(f,a,b,n)
   exactm = tril(exact*ones(n,n));
   relerr = (tableau-exactm)/exact
end


% Sauer's Program 5.1 Romberg integration
% Computes approximation to definite integral
% Inputs: Matlab inline function specifying integrand f,
%    a,b integration interval, n=number of rows
% Output: Romberg tableau r
function r=romberg(f,a,b,n)
  h=(b-a)./(2.^(0:n-1));
  r(1,1)=(b-a)*(f(a)+f(b))/2;
  for j=2:n
    subtotal = 0;
     for i=1:2^(j-2)
      subtotal = subtotal + f(a+(2*i-1)*h(j));
   end
    r(j,1) = r(j-1,1)/2+h(j)*subtotal;
    for k=2:j
       r(j,k) = (4^(k-1)*r(j,k-1)-r(j-1,k-1))/(4^(k-1)-1);
    end
  end
end

>> useromberg
n =
     6
f = 
    @exp
a =
     0
b =
     1
tableau =
  Columns 1 through 5 
     1.859140914229523e+00                         0                         0                         0                         0
     1.753931092464825e+00     1.718861151876593e+00                         0                         0                         0
     1.727221904557517e+00     1.718318841921747e+00     1.718282687924758e+00                         0                         0
     1.720518592164302e+00     1.718284154699897e+00     1.718281842218440e+00     1.718281828794531e+00                         0
     1.718841128579994e+00     1.718281974051892e+00     1.718281828675358e+00     1.718281828460389e+00     1.718281828459078e+00
     1.718421660316328e+00     1.718281837561772e+00     1.718281828462431e+00     1.718281828459051e+00     1.718281828459046e+00
  Column 6 
                         0
                         0
                         0
                         0
                         0
     1.718281828459046e+00
exact =
     1.718281828459046e+00
relerr =
  Columns 1 through 5 
     8.197670686932633e-02                         0                         0                         0                         0
     2.074704126839902e-02     3.371527347564584e-04                         0                         0                         0
     5.202916046949473e-03     2.154097313295723e-05     5.001890248271996e-07                         0                         0
     1.301744375229878e-03     1.353817990013009e-06     8.007647236210882e-09     1.952444398693568e-10                         0
     3.254996425414904e-04     8.473164531811413e-08     1.258889792761678e-10     7.816807423238881e-13     1.899604382899844e-14
     8.137888381643365e-05     5.297574834199787e-09     1.970161117121838e-12     3.230619698809258e-15     2.584495759047406e-16
  Column 6 
                         0
                         0
                         0
                         0
                         0
     2.584495759047406e-16
>> 

>> useromberg
n =
     8
f = 
    @(x) 1/(1+12*x^2)
a =
     0
b =
     1
exact =
   3.7232e-01
tableau =
   5.3846e-01            0            0            0            0            0            0            0
   3.9423e-01   3.4615e-01            0            0            0            0            0            0
   3.7223e-01   3.6490e-01   3.6615e-01            0            0            0            0            0
   3.7214e-01   3.7211e-01   3.7259e-01   3.7269e-01            0            0            0            0
   3.7228e-01   3.7232e-01   3.7234e-01   3.7233e-01   3.7233e-01            0            0            0
   3.7231e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01            0            0
   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01            0
   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01
exactm =
   3.7232e-01            0            0            0            0            0            0            0
   3.7232e-01   3.7232e-01            0            0            0            0            0            0
   3.7232e-01   3.7232e-01   3.7232e-01            0            0            0            0            0
   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01            0            0            0            0
   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01            0            0            0
   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01            0            0
   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01            0
   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01   3.7232e-01
relerr =
   4.4623e-01            0            0            0            0            0            0            0
   5.8843e-02  -7.0284e-02            0            0            0            0            0            0
  -2.4565e-04  -1.9942e-02  -1.6586e-02            0            0            0            0            0
  -4.9422e-04  -5.7708e-04   7.1391e-04   9.8851e-04            0            0            0            0
  -1.2409e-04  -7.0549e-07   3.7720e-05   2.6986e-05   2.3216e-05            0            0            0
  -3.1035e-05  -1.8897e-08   2.6876e-08  -5.7142e-07  -6.7949e-07  -7.0285e-07            0            0
  -7.7598e-06  -1.1831e-09  -2.2166e-12  -4.2886e-10   1.8103e-09   2.4763e-09   2.6486e-09            0
  -1.9400e-06  -7.3977e-11  -3.4888e-14  -2.9819e-16   1.6815e-12  -8.6624e-14  -6.9135e-13  -8.5312e-13
>>