function [q,cnt] = romberg(a, b, tol)

%  Approximates the integral from a to b of f(x)dx to absolute tolerance
%  tol by using the trapezoidal rule with repeated Richardson extrapolation.
%  Returns number of function evaluations in cnt.

%  Make first estimate using one interval.

n = 1;  h = b-a;
fv = [f(a); f(b)];
est(1,1) = .5*h*(fv(1)+fv(2));
cnt = 2;                         % Count no. of function evaluations.

%  Keep doubling the number of subintervals until desired tolerance is
%  achieved or max no. of subintervals (2^10 = 1024) is reached.

err = 100*abs(tol);              % Initialize err to something > tol.

k = 0;
while err > tol & k < 10,
  k = k+1;  n = 2*n;  h = h/2;
  fvnew = zeros(n+1,1);          % Store computed values of f to reuse
  for i=1:2:n+1,                 % when h is cut in half.
    fvnew(i) = fv((i-1)/2 + 1);
  end;
  for i=2:2:n,                   % Compute f at midpoints of previous intervals
    fvnew(i) = f(a+(i-1)*h);
  end;
  cnt = cnt + (n/2);             % Update no. of function evaluations.
  fv = fvnew;
  trap = .5*(fv(1)+fv(n+1));     % Use trapezoidal rule with new h value
  for i=2:n,                     % to estimate integral.
    trap = trap + fv(i);
  end;
  est(k+1,1) = h*trap;         % Store new estimate in first column of tableau.

%    Perform Richardson extrapolations.

  for j=2:k+1,
    est(k+1,j) = ((4^(j-1))*est(k+1,j-1) - est(k,j-1))/(4^(j-1)-1);
  end;
  q = est(k+1,k+1);

%    Estimate error.  (This is usually an overestimate of the error in q.
%    It is a better estimate of the error at previous stage.)
  if k > 2, err = abs(q - est(k,k)); end;

%    Print out approximation to integral and error at each step,
%    for monitoring convergence.
  cnt, q, err

end;