/*     Newton's Method to Solve a Nonlinear Function      */
   /* The user must supply initial values,                   */
   /* and the FUN and DERIV functions.                       */
   /* On entry: FUN evaluates the function f in terms of x   */
   /* initial values are given to x                          */
   /* DERIV evaluates jacobian j                             */
   /* tuning variables: CONVERGE, MAXITER.                   */
   /* On exit: solution in x, function value in f close to 0 */
   /* ITER has number of iterations.                         */
*ods trace output;
proc iml;
start newton;
   run fun;          /* evaluate function at starting values */
   do iter = 1 to maxiter           /* iterate until maxiter */
   while(max(abs(f))>converge); /* iterations or convergence */
      run deriv;                /* evaluate derivatives in j */
      delta = -solve(j,f);    /* solve for correction vector */
      x = x+delta;                  /* the new approximation */
      run fun;                      /* evaluate the function */
   end;
finish newton;

maxiter = 15;                  /* default maximum iterations */
converge = .000001;         /* default convergence criterion */

   /* User-supplied function evaluation */
start fun;
   x1 = x[1] ;
   x2 = x[2] ;                         /* extract the values */
   f = (x1+x2-x1*x2+2) //
       (x1*exp(-x2)-1);             /* evaluate the function */
finish fun;

   /* User-supplied derivatives of the function */
start deriv;
   /* evaluate jacobian */
   j = ((1-x2)||(1-x1) ) // (exp(-x2)||(-x1*exp(-x2)));
finish deriv;

do;
   print "Solving the system: X1+X2-X1*X2+2=0, X1*EXP(-X2)-1=0" ,;
   x={.1, -2};     /* starting values */
   run newton;
   print x f;
end;