SUBROUTINE NEWTON(F, DF, XINIT, EPS, ITMAX, ROOT, IER)
C     -----------------
C
C     THIS ROUTINE CALCULATES A ROOT OF
C               F(X) = 0
C     USING NEWTON'S METHOD.
C
C     INPUT PARAMETERS:
C     F,DF:  THE NAMES OF FUNCTION SUBPROGRAMS FOR COMPUTING
C            F(X) AND DF(X)=F'(X).  THE ACTUAL NAMES USED IN
C            CALLING SUBROUTINE NEWTON MUST BE DECLARED IN
C            AN EXTERNAL STATEMENT IN THE CALLING PROGRAM.
C     XINIT: AN INITIAL GUESS OF THE SOLUTION.
C     EPS:   THE DESIRED ERROR TOLERANCE.  THE TEST FOR
C            CONVERGENCE IS
C                    ABS(X1-X0) .LE. EPS,
C            WHERE X1 IS THE CURRENT NEWTON ITERATE AND X0
C            IS THE PRECEDING ITERATE.
C     ITMAX: AN UPPER LIMIT ON THE NUMBER OF ITERATES
C            TO BE COMPUTED. WHEN THE NUMBER OF ITERATES
C            REACHES THIS NUMBER, THE PROGRAM IS TERMINATED.
C            WARNING: ITMAX IS ALTERED ON OUTPUT.
C
C     OUTPUT PARAMETERS:
C     ROOT:  CONTAINS THE COMPUTED VALUE OF THE SOLUTION.
C            REGARDLESS OF HOW THE ROUTINE IS TERMINATED,
C            SUCCESSFULLY OR OTHERWISE, ROOT WILL CONTAIN
C            THE MOST RECENTLY COMPUTED NEWTON ITERATE.
C     ITMAX: SET TO THE NUMBER OF ITERATES COMPUTED IN
C            THE ROUTINE.
C     IER:   AN ERROR INDICATOR.
C        =0  MEANS A SUCCESSFUL COMPLETION OF NEWTON.
C        =1  MEANS ITMAX ITERATES WERE COMPUTED AND THE
C            ROUTINE WAS ABORTED.
C        =2  MEANS THAT THE DERIVATIVE DF(X) BECAME ZERO
C            AT SOME NEWTON ITERATE, AND THE ROUTINE
C            WAS ABORTED.
C
      PARAMETER(ZERO=0.0)
C
C     INITIALIZE.
      X0 = XINIT
      ITNUM = 1
C
C     BEGIN MAIN LOOP.
10    DENOM = DF(X0)
      IF(DENOM .EQ. ZERO) THEN
C         DERIVATIVE EQUALS ZERO. TERMINATE ITERATION.
          IER = 2
          ROOT = X0
          ITMAX = ITNUM-1
          RETURN
      END IF 
C
C     COMPUTE NEWTON ITERATION.
      X1 = X0 - F(X0)/DENOM 
C
      IF(ABS(X1-X0) .LE. EPS) THEN
C         ERROR TEST SATISFIED.
          IER = 0
          ROOT = X1
          ITMAX = ITNUM
          RETURN
      END IF
C
      IF(ITNUM .LE. ITMAX) THEN
C         INITIALIZE FOR ANOTHER LOOP.
          ITNUM = ITNUM + 1
          X0 = X1
          GO TO 10
      ELSE
C         ITMAX ITERATES HAVE BEEN COMPUTED. TERMINATE.
          IER = 1
          ROOT = X1
          RETURN
      END IF
      END