/****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: NEWTON                                              */
 /*   TITLE: NEWTON'S METHOD FOR SOLVING NONLINEAR SYSTEMS       */
 /* PRODUCT: IML                                                 */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS: MATRIX  SUGI6                                       */
 /*   PROCS: IML                                                 */
 /*    DATA:                                                     */
 /*                                                              */
 /* SUPPORT: RHD                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/


proc iml;


 /*------ 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 zero  */
 /*           ITER has number of iterations.                    */
 /*-------------------------------------------------------------*/
START NEWTON;
  RUN FUN;                /* EVALUATE FUNCTION AT STARTING VALUES */
  DO ITER=1 TO MAXITER    /* ITERATE UNTIL MAXITER ITERATIONS */
     WHILE(MAX(ABS(F))>CONVERGE); /* 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;

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);     /* EVAULATE THE FUNCTION  */
FINISH;


 /*---user-supplied derivatives of the function---*/
START DERIV;
  J=(   (1-X2)||(1-X1)        )//
    ( EXP(-X2)||(-X1*EXP(-X2)));  /* EVALUATE JACOBIAN */
FINISH;

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;