newton := proc(  f :: algebraic,
                 x :: name,
                x0 :: numeric,
               tol :: numeric
                )
   local g, z, k;

   g := diff(f, x);  # get the derivative
   z := Vector(20);  # open an empty vector of dimension 20
   z[1] := x0;       # initialize the initial iterate

   for k from 2 to 20 do

      z[k] := evalf(z[k-1] - subs(x=z[k-1],f)/subs(x=z[k-1],g));

      print(k,z[k],abs(z[k]-z[k-1]));

      if abs(z[k]-z[k-1]) < tol then
         return z[k];
      end if;

   end do;

   print("Newton iteration failed");

end proc;