      PROGRAM BISECT
************************************************************************
* Program to find an approximate zero of a function F in an interval   *
* [A, B} using the bisection method.  Identifiers used:                *
*     A, B   : [A,B] is the interval in which a root is to be found    *
*     MID    : midpoint of subinterval being searched for a root       *
*     FA, FB : value of function F at A and B, respectively            *
*     FMID   :   "    "    "     "  " MID                              *
*     ZERO   : true if exact zero found, else false                    *
*     INTLEN : length of current subinterval                           *
*     EPSIL  : terminate when INTLEN < EPSIL                           *
*     F      : the function whose root is being sought                 *
*                                                                      *
* Input:  Endpoints A and B and desired accuracy EPSIL                 *
* Output: Root or approximate root or message indicating that interval *
*         [A, B} may not contain a root                                *
************************************************************************

      REAL F, A, B, MID, FA, FB, FMID, INTLEN, EPSIL
      LOGICAL ZERO

      PRINT *, 'ENTER INTERVAL ENDPOINTS A AND B, AND EPSILON:'
      READ *, A, B, EPSIL
      FA = F(A)
      FB = F(B)
      IF (FA .EQ. 0) THEN
         PRINT *, 'ROOT = ', A
      ELSE IF (FB .EQ. 0) THEN
         PRINT *, 'ROOT = ', B
      ELSE IF (FA * FB .GT. 0) THEN
         PRINT *, 'F(A) AND F(B) HAVE SAME SIGN.'
         PRINT *, 'INTERVAL MAY NOT CONTAIN A ROOT.'
      ELSE
         ZERO = .FALSE.
         INTLEN = B - A

*        While no exact zero found and length of interval is greater
*        than EPSIL, bisect the interval and search one of the halves
10       IF (.NOT. ZERO  .AND. (INTLEN / 2.0 .GT. EPSIL)) THEN
            MID = (A + B) / 2.0
            FMID = F(MID)
            IF (FMID .EQ. 0) THEN
               ZERO = .TRUE.
            ELSE IF (FA * FMID .LT. 0) THEN
*              Root in left half -- search it
               B = MID
               FB = FMID
            ELSE

*              Root in right half -- search it
               A = MID
               FA = FMID
            END IF
            INTLEN = INTLEN/2
         GO TO 10
         END IF
         IF (ZERO) THEN
            PRINT *, 'ROOT =', MID
         ELSE
            PRINT *, 'APPROX. ROOT =   ', MID
            PRINT *, 'FUNCTION VALUE = ', FMID
         END IF
      END IF
      END

** F(X) *****************************************
* The function for which a root is being sought *
*************************************************

      REAL FUNCTION F(X)
      REAL X

      F = X**3 + X - 5.0

      END

