function x = newton_root (f, df, x0, tol, maxiter)

  % f:        a given scalar function y=f(x)
  % df:       the derivative of f
  % x0:       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 = x0;
  xold = inf;
  j = 0;
  while (abs(xold-xnew) > tol &
	 j < maxiter)
    j = j+1;
    xold = xnew;
    xnew = xold - feval(f,xold) / feval(df,xold);
  end
  x = xnew;

end

  