function x = secant_root (f, x0, x1, tol, maxiter)

  % f:        a given scalar function y=f(x)
  % x0:       first starting point
  % x1:       second starting point
  % tol:      stop when the difference of two iterates is less than tol
  % maxiter:  stop after at most maxiter iterations in any case

  xnew = x1;
  xold = x0;
  j = 0;
  while (abs(xold-xnew) > tol &
	 j < maxiter)
    j = j+1;
    xtmp = xnew;
    xnew = xnew - feval(f,xnew) * ...
                  (xold - xnew) / (feval(f,xold) - feval(f,xnew));
    xold = xtmp;
  end
  x = xnew;

end

  