x={.2,.2,.1,.1,.3};
u=x[1]^.3*x[2]^.2*x[3]^.1*x[4]^.2*x[5]^.2;
a={10 3 2 7 5};
b={100};

proc g(x);
   local du,u;
   u=x[1]^.3*x[2]^.2*x[3]^.1*x[4]^.2*x[5]^.2;
   du=u*(.3/x[1]|.2/x[2]|.1/x[3]|.2/x[4]|.2/x[5]);
   retp(du);
   endp;

proc h(x);
   local ddu,u;
   u=x[1]^.3*x[2]^.2*x[3]^.1*x[4]^.2*x[5]^.2;
   ddu=u*( (.3*(.3-1))/(x[1]*x[1])~   (.3*.2)/(x[1]*x[2]) ~  (.3*.1)/(x[1]*x[3])  ~   (.3*.2)/(x[1]*x[4])  ~  (.3*.2)/(x[1]*x[5])|
             (.3*.2)/(x[2]*x[1])  ~(.2*(.2-1))/(x[2]*x[2])~  (.2*.1)/(x[2]*x[3])  ~   (.2*.2)/(x[2]*x[4])  ~  (.2*.2)/(x[2]*x[5])|
             (.3*.1)/(x[3]*x[1])  ~   (.2*.1)/(x[3]*x[2]) ~(.1*(.1-1))/(x[3]*x[3])~   (.1*.2)/(x[3]*x[4])  ~  (.1*.2)/(x[3]*x[5])|
             (.3*.2)/(x[4]*x[1])  ~   (.2*.2)/(x[4]*x[2]) ~  (.2*.1)/(x[4]*x[3])  ~ (.2*(.2-1))/(x[4]*x[4])~  (.2*.2)/(x[4]*x[5])|
             (.3*.2)/(x[5]*x[1])  ~   (.2*.2)/(x[5]*x[2]) ~  (.2*.1)/(x[5]*x[3])  ~   (.2*.2)/(x[5]*x[4])  ~(.2*(.2-1))/(x[5]*x[5]));
   retp(ddu);
   endp;

proc step(p,x);
   local alpha,s;
   s=((b-a*x)./(a*p));
   if s >= 1;
      alpha=1;
   else;
      alpha=s;
   endif;
   retp(alpha);
   endp;

icount=1;
do until icount > 20;
   ei=eigrs(h(x));
   if maxc(ei) < -0.001;
      hh=inv(h(x));
      print "Newton";
   else;
      hh=eye(5);
   endif;
   p=inv(hh)*g(x);
   s=step(p,x);
   if s == 1;
      x=x+s*p;
      icount=icount+1;
   else;
      x=x+s*p;
      icount=21;
   endif;
   u=x[1]^.3*x[2]^.2*x[3]^.1*x[4]^.2*x[5]^.2;
   print u~icount~s~p';
   endo;

z=null(a);
icount=1;
do until icount > 20;
   hh=z'h(x)*z;
   ei=eigrs(hh);
   if maxc(ei) < -0.001;
      p=-inv(hh)*(z'g(x));
      x=x+z*p;
      ty=1;
   else;
      x=x+z*z'g(x);
      ty=-1;
   endif;
   u=x[1]^.3*x[2]^.2*x[3]^.1*x[4]^.2*x[5]^.2;
   u~icount~(z'g(x))'~ty;
   icount=icount+1;
   test=sqrt((z'g(x))'(z'g(x)));
   if test < 0.00001;
      icount=1000;
   endif;
   endo;