      PROGRAM NEWTON
************************************************************************
* Program to find an approximate root of a function F using Newton's   *
* method.  Variables used are:                                         *
*     XOLD   : previous approximation (initially the first one)        *
*     FPOLD  : value of the derivative of F at XOLD                    *
*     XNEW   : the new approximation                                   *
*     FVALUE : value of F at an approximation                          *
*     EPSIL  : repetition stops when ABS(FVALUE) is less than EPSIL    *
*     NUMITS : limit on number of iterations                           *
*     N      : number of iterations                                    *
*                                                                      *
* Input:  EPSIL, NUMITS, and XOLD                                      *
* Output: A sequence of approximations to the root or an error message *
*         indicating that the method fails                             *
************************************************************************

      INTEGER NUMITS, N
      REAL XOLD, FPOLD, XNEW, EPSIL, F, FPRIME, FVALUE

* Get termination values EPSIL & NUMITS and initial approximation

      PRINT *, 'ENTER EPSILON, LIMIT ON # OF ITERATIONS, AND'
      PRINT *, 'THE INITIAL APPROXIMATION'
      READ *, EPSIL, NUMITS, XOLD

* Initialize FVALUE and N; print headings and initial values

      FVALUE = F(XOLD)
      N = 0
      PRINT *, '  N      X(N)        F(X(N))'
      PRINT *, '============================='
      PRINT 10, 0, XOLD, FVALUE
 10   FORMAT (1X, I3, F11.5, E14.5)

* Iterate using Newton's method while ABS(FVALUE) is greater
* than or equal to EPSIL and N has not reached NUMITS
* Terminate if the derivative is 0 at some approximation

 20   IF ((ABS(FVALUE) .GE. EPSIL) .AND. (N .LT. NUMITS)) THEN
         N = N + 1
         FPOLD = FPRIME(XOLD)
         IF (FPOLD .EQ. 0) THEN
            PRINT *, 'NEWTON''S METHOD FAILS -- DERIVATIVE = 0'
*           Force repetition to terminate
            N = NUMITS
         ELSE
            XNEW = XOLD - (FVALUE / FPOLD)
            FVALUE = F(XNEW)
            PRINT 10, N, XNEW, FVALUE
            XOLD = XNEW
         END IF
      GO TO 20
      END IF

      END

**F(X)**************************************
* Function for which a root is being found *
********************************************

      FUNCTION F(X)

      REAL X, F

      F = X**3 + X - 5
      END

**FPRIME(X)*************************
* The derivative of the function f *
************************************

      FUNCTION FPRIME(X)

      REAL X, FPRIME

      FPRIME = 3*X**2 + 1
      END
