function r = secant( fun, p0, p1, tol, nmax, verbose)

% secant  Uses the secant method to find a root 
%         of the equation  fun(x) = 0 
%         starting at with the points p0 and p1
%
% Input:  fun = the function whose zero we are searching
%         p0, p1 = initial approximations
%         tol = x tolerance
%         nmax = maximum number of iterations
%         verbose = 1 in case all iterations are needed, 
%                   0 otherwise
%
% Output: r = the root of fun(x) = 0 with the desired tolerance

if tol < 3*eps       %  The tolerance is too small
  warning(sprintf('the tolerance is too small\n'));
  return;
end

if verbose == 1         %  Printing the headings
  fprintf(['\niters         x                 x-xold' ...
	   '    log_10(|x-xold|)\n\n']);
end


%  The main part of the program starts here

pold = p0;
pnew = p1;
fold = feval(fun,pold);
fnew = feval(fun,pnew);

for i = 1:nmax
  pnext = pnew - fnew*(pnew-pold)/(fnew-fold);  

  if verbose == 1             %  Printing all iterations
    fprintf('%3d  %18.14f  %18.14f  %10.5f\n', ...
		        i, pnext, pnext-pnew,log(abs(pnext-pnew))/log(10.) );
  end

  if abs(pnext-pold) < tol
    r = pnext;    return;      %  We found the root
  end

  pold = pnew;
  fold = feval(fun,pold);
  pnew = pnext;
  fnew = feval(fun,pnew); 

end

warning(sprintf('root not found within tol after %d steps\n',i));