function [I,T]=Romberg(f,a,b,tol);
% ROMBERG quadrature using the Romberg scheme 
%   [I,T]=Romberg(f,a,b,tol) computes an approximation of int_a^b f(x)
%   dx to a relative tolerance tol using the Romberg scheme.

h=b-a; intv=1; s=(f(a)+f(b))/2; 
T(1,1)=s*h;   
for i=2:15
  intv=2*intv; h=h/2;
  s=s+sum(f(a+[1:2:intv]*h));
  T(i,1)=s*h; 
  vhj=1;
  for j=2:i
    vhj=vhj/4; 
    T(i,j)=(vhj*T(i-1,j-1)-T(i,j-1))/(vhj-1);
  end;
  if abs(T(i,i)-T(i-1,i-1))<tol*abs(T(i,i)),
    I=T(i,i); return
  end
end 
warning(['limit of extrapolation steps reached. ',...
         'Required tolerance may not be met.']);
I=T(i,i);