function z = buelerromberg(f,a,b,n)
% BUELERROMBERG  Approximate integral
%    /b
%    |   f(x) dx
%    /a
% using Romberg integration with any number (n) of levels.
% The function f(x) must be vectorized (i.e. an input list x must be allowed).
%
% Example:  >> format long
%           >> buelerromberg(@(x) exp(-x),0,1,3)
%           >> exact = 1 - exp(-1)
%           >> buelerromberg(@(x) cos(x).^3,-1,2,7)
%           >> F = @(u) u - u^3/3; exact = F(sin(2)) - F(sin(-1))

if n < 1, error('n must be 1,2,3,...'), end

h = b - a;                              % initial h
R = zeros(n,n);                         % space for whole table
R(1,1) = 0.5 * h * (feval(f,a) + feval(f,b));  % n = 1 version of trapezoid
L = 1;                                  % L = 2^{i-2}
for i=2:n
  s = 0.0;
  for k = 1:L
    s = s + feval(f,a + (k-0.5)*h);
  end
  R(i,1) = 0.5 * (R(i-1,1) + h * s);    % n = 2^{i-1} case of trapezoid
  for j=2:i
    c = 1.0 / (4^(j-1) - 1);
    R(i,j) = R(i,j-1) + c * ( R(i,j-1) - R(i-1,j-1) );
  end
  h = h / 2;                            % h value for next iteration
  L = L * 2;                            % L value for next iteration
end
z = R(n,n);
% R                                     % uncomment to see whole table R_{k,j}
