# 
# Program of Newton's iteration for
#   solving the equation f(x)=0
#
#  Input:    f --- the function 
#                    (e.g.  [>f:=x-> x^2-2;  )
#           x0 --- the initial iterate
#            n --- the maximum number of iterative 
#                     steps allowed
#          tol --- the error tolerance
#
newton:=proc(f,x0,n,tol)
   local x, k, s, g, delta;

   x:=array(0..n);     # open an empty array for iterates
   g:=D(f);            # differentiate the function f
   x[0]:=x0;           # initialize the 0-th iterate
   delta:=tol+1;       # initialize the error estimate
                       
                       # the loop of Newton's iteration
   for k to n while abs(delta)>tol do
      delta:= evalf( f(x[k-1])/g(x[k-1]) ); 
      s:=x[k-1]-delta;
      x[k]:=s;
   od;
                       # checking convergence
   if abs(delta)<=tol then
      print(`The solution`);
      print(s);
   else
      print(`no convergence yet`);
   fi;

end;