
      REAL FUNCTION AFunction (X)
      REAL X
        AFunction = (2 * X + 3) * (X - 3)
      RETURN
      END !FUNCTION

      REAL FUNCTION AFunctionPrime (X)
      REAL X
        AFunctionPrime = 4 * X - 3
      RETURN
      END  !FUNCTION




      FUNCTION BisectionRoots (x1, x2, Tolerance, maxIterations,
     +        valueAtRoot, errr)
      REAL x1, x2, Tolerance, valueAtRoot, FPP, nearzero
      INTEGER maxIterations, errr

      PARAMETER (FPP = 1E-15)
      PARAMETER (nearzero = 1E-20 )
      REAL result, FA,  FC, XA, XB, XC
      INTEGER i
      LOGICAL done


      errr = 0
      i = 0
      done = .FALSE.
      XA = x1
      XB = x2
      FA = AFunction(XA)
      FB = AFunction(XB)
      IF ( RootBracketed(FA, FB) .EQ. 0 ) THEN
        errr = 1
      ELSE
        DO WHILE (.NOT. done .AND. i .NE. maxIterations)
          XC = (XA + XB) / 2.0
          FC = AFunction(XC)
          IF (RootBracketed(FA, FC) .EQ. 1) THEN
            XB = XC
          ELSE
            XA = XC
            FA = FC
          END IF
          IF (ABS(XB - XA) .LE. Tolerance) THEN
            done = .TRUE.
            result = (XA + XB) / 2.0
            valueAtRoot = AFunction(result)
          END IF
          i = i + 1
        END DO
        IF (i .EQ. maxIterations) errr = 2
      END IF
      BisectionRoots = result
      END !FUNCTION



      REAL FUNCTION BrentRoots (x1, x2, Tolerance, maxIterations,
     +    valueAtRoot, errr)
      REAL x1, x2, Tolerance, valueAtRoot, FPP, nearzero
      INTEGER maxIterations, errr

      PARAMETER (FPP = 1E-15)
      PARAMETER (nearzero = 1E-20)
      REAL result, AA, BB, CC, DD, EE, FA, FB
      REAL  Tol1, PP, QQ, RR, SS, XM
      INTEGER i
      LOGICAL done
      i = 0
      done = .FALSE.
      errr = 0
      AA = x1
      BB = x2
      FA = AFunction(AA)
      FB = AFunction(BB)
      IF ( RootBracketed(FA, FB) .EQ. 0) THEN
         errr = 1
      ELSE
        FC = FB
        DO WHILE (.NOT. done .AND. i .NE. maxIterations)
          IF  (RootBracketed(FC, FB) .EQ. 0) THEN
            CC = AA
            FC = FA
            DD = BB - AA
            EE = DD
          END IF
          IF (ABS(FC) .LT. ABS(FB)) THEN
            AA = BB
            BB = CC
            CC = AA
            FA = FB
            FB = FC
            FC = FA
          END IF
          Tol1 = 2.0 * FPP * ABS(BB) + 0.5 * Tolerance
          XM = 0.5 * (CC - BB)
          IF ((ABS(XM) .LE. Tol1) .OR. (ABS(FA) .LT. nearzero)) THEN
            result = BB
            done = .TRUE.
            valueAtRoot = AFunction(result)
          ELSE
            IF ((ABS(EE) .GE. Tol1) .AND. (ABS(FA) .GT. ABS(FB))) THEN
              SS = FB / FA
              IF (ABS(AA - CC) .LT. nearzero) THEN
                PP = 2.0 * XM * SS
                QQ = 1.0 - SS
              ELSE
                QQ = FA / FC
                RR = FB / FC
                PP = SS * (2.0 * XM * QQ *
     +               (QQ - RR) - (BB - AA) * (RR - 1.0))
                QQ = (QQ - 1.0) * (RR - 1.0) * (SS - 1.0)
              END IF
              IF (PP .GT. nearzero)  QQ = -QQ
              PP = ABS(PP)
              IF (2.0 * PP .LT.
     +           Minimum(3.0*XM*QQ - ABS(Tol1*QQ), ABS(EE*QQ))) THEN
                EE = DD
                DD = PP / QQ
              ELSE
                DD = XM
                EE = DD
              END IF
            ELSE
              DD = XM
              EE = DD
            END IF
            AA = BB
            FA = FB
            IF (ABS(DD) .GT. Tol1) THEN
              BB = BB + DD
            ELSE
              IF (XM .GT. 0) THEN
                BB = BB + ABS(Tol1)
              ELSE
                BB = BB - ABS(Tol1)
              END IF
            END IF
            FB = AFunction(BB)
          END IF
          i = i + 1
        END DO
        IF (i .EQ. maxIterations)  errr = 2
      END IF
      BrentRoots = result
      END !FUNCTION



      FUNCTION Minimum (x1, x2)
      REAL x1, x2
      REAL result

      IF (x1 .LT. x2) THEN
        result = x1
      ELSE
        result = x2
      END IF
      Minimum = result
      END !FUNCTION



      FUNCTION NewtonRoots (x1, Tolerance, maxIterations,
     +     valueAtRoot, errr)
      REAL FPP, nearzero
      PARAMETER ( FPP = 1E-15)
      PARAMETER ( nearzero = 1E-20)
      REAL result, lastX, FX, FPX, X
      LOGICAL done
      INTEGER maxIterations, errr

      i = 0
      errr = 0
      done = .FALSE.
      X = x1
      FX = AFunction(X)
      DO WHILE (.NOT. done .AND. i .LE. maxIterations .AND. errr .EQ. 0)

       IF (ABS(X) .LT. nearzero) THEN
         errr = 1
       ELSE
         FPX = AFunctionPrime(X)
         IF (ABS(FPX) .LT. nearzero) THEN
           errr = 2
         ELSE
           lastX = X
           X = X - FX / FPX
           FX = AFunction(X)
           IF (ABS(X - lastX) .LE. Tolerance) THEN
             done = .TRUE.
             result = X
             valueAtRoot = AFunction(X)
           END IF
           i = i + 1
         END IF
        END IF
      END DO !! UNTIL done .OR. (i .GT. maxIterations) .OR. (errr .NE. 0)
      IF( i .EQ. maxIterations)  errr = 3
      NewtonRoots = result

      END ! FUNCTION


      FUNCTION RootBracketed (x1, x2)
      REAL x1, x2
      REAL result
      IF ( (x1 .GT. 0.0 .AND. x2 .GT. 0.0) .OR.
     +   (x1 .LT. 0.0 .AND. x2 .LT. 0.0)) THEN
        result = 0
      ELSE
        result =  1
      END IF
      RootBracketed = result
      END !FUNCTION

