** SECANT **************************************************************
* Subroutine to find an approximate root of a function F using the     *
* secant method.  The argumnets are:                                   *
*     X0, X1  : two initial approximations                             *
*     NUMITS  : limit on the number of iterations                      *
*     EPSIL   : repetition stops when ABS(FVALUE) is less than EPSIL   *
* Local variables are:                                                 *
*     N       : number of iterations                                   *
*     XNEW    : the new approximation                                  *
*     F0, F1  : values of F at X0 and X1, respectively                 *
*     FNEW    : value of F at XNEW                                     *
*                                                                      *
* Accepts: X0, X1, NUMITS, and EPSIL                                   *
* Output:  A sequence of approximations to the root or an error        *
*          message indicating that the method fails                    *
************************************************************************

      SUBROUTINE SECANT(X0, X1, NUMITS, EPSIL)

      INTEGER NUMITS, N
      REAL X0, X1, EPSIL, F, XNEW, F0, F1, FNEW

      PRINT *, '  N      X(N)    ABS(F(X(N)))'
      PRINT *, '============================='
      F0 = F(X0)
      F1 = F(X1)
      PRINT *, 0, X0, ABS(F0)
      PRINT *, 1, X1, ABS(F1)

* Iterate using the secant method while ABS(FNEW) is greater
* than or equal to EPSIL and N has not reached NUMITS.
* Terminate if F0 = F1.

      N = 1
10    IF ((ABS(F1) .GE. EPSIL) .AND. (N .LT. NUMITS)) THEN
          N = N + 1
          IF (F0 .EQ. F1) THEN
              PRINT *, 'SECANT METHOD FAILS'
*             Force repetition to terminate
              N = NUMITS
          ELSE
              XNEW = X1 - F(X1) * ((X1 - X0) / (F1 - F0))
              FNEW = F(XNEW)
              PRINT *, N, XNEW, ABS(FNEW)
              X0 = X1
              F0 = F1
              X1 = XNEW
              F1 = FNEW
          END IF
      GO TO 10
      END IF
      END
