NewtonRaphson := proc(f,x0,tol,maxiter)
  
  #
  # Newton-Raphson iteration to find a root of f
  #
  # usage:  NewtonRaphson(f,x0,tol,maxiter)
  # where:  f = the function in the equation f(x) = 0
  #         x0 = the starting point
  #         tol = tolerance (iteration continues until
  #               two steps give values differing 
  #               by less than that amount)
  #         maxiter = maximum number of steps that
  #                   will be performed
  #
  
    local x,xnext,i;
    
    x:=x0;
    for i to maxiter do
       xnext := evalf(x - f(x)/D(f)(x));
       if abs(xnext - x) < tol then
          return xnext;
       end if;
       x := xnext; 
    end do;
    lprint("Newton-Raphson iteration did not converge in",maxiter,"steps");
    end proc: