function [val,err,ampl] = richardson(tol,p,nmin,f,varargin)

% [val,err,ampl] = richardson(tol,p,nmin,f,varargin)
%
% richardson extrapolation with error control
%
%   tol     requested relative error
%   p       f(...,h) has an asymptotic expansion in h^p
%   nmin    begin double harmonic sequence at n=2*nmin
%   f       function, evaluates as f(varargin,h)
%
%   val     extrapolated value (h -> 0)
%   err     estimated relative error
%   ampl    amplification of relative error in the 
%               evaluation of f

% initialize tableaux
j_max = 8; T = zeros(j_max,j_max);
n = 2*(nmin:(j_max+nmin+1));        % double harmonic sequence
j=1; h=1/n(1);
T(1,1) = feval(f,varargin{:},h);
val=T(1,1); err=abs(val); 

% do Richardson until convergence
while err(end)/abs(val) > tol & j<j_max
    j=j+1; h=1/n(j);
    T(j,1) = feval(f,varargin{:},h);
    A(j,1) = abs(T(j,1));
    for k=2:j
        T(j,k) = T(j,k-1) + (T(j,k-1)-T(j-1,k-1))/((n(j)/n(j-k+1))^p-1);
        A(j,k) = A(j,k-1) + (A(j,k-1)+A(j-1,k-1))/((n(j)/n(j-k+1))^p-1);
    end
    val = T(j,j);
    err = abs(diff(T(j,1:j)));   % subdiagonal error estimates
    if j > 2            % extrapolate error estimate once more
        err(end+1) = err(end)^2/err(end-1);
    end
end
ampl = A(end,end)/abs(val);
err = max(err(end)/abs(val),ampl*eps);

return;