      FUNCTION ZEROIN(F,AX,BX,AERR,RERR)
      REAL F,AX,BX,AERR,RERR
      EXTERNAL F
C
C     A ZERO OF THE FUNCTION F(X) IS COMPUTED IN THE INTERVAL (AX,BX).
C
C  INPUT..
C
C  F      FUNCTION SUBPROGRAM WHICH EVALUATES F(X) FOR ANY X IN
C         THE CLOSED INTERVAL (AX,BX)
C  AX     LEFT ENDPOINT OF THE INTERVAL
C  BX     RIGHT ENDPOINT OF THE INTERVAL
C  AERR   THE ABSOLUTE ERROR TOLERANCE TO BE SATISFIED
C  RERR   THE RELATIVE ERROR TOLERANCE TO BE SATISFIED
C
C
C  OUTPUT..
C
C  ZEROIN ABCISSA APPROXIMATING A ZERO OF F IN THE INTERVAL (AX,BX)
C
C
C      IT IS ASSUMED  THAT  F(AX)  AND  F(BX)  HAVE  OPPOSITE  SIGNS
C  WITHOUT  A  CHECK.  ZEROIN  RETURNS A ZERO X IN THE GIVEN INTERVAL
C  (AX,BX) TO WITHIN A TOLERANCE  MAX(RERR,4*EPS)*ABS(X) + AERR, WHERE
C  EPS IS THE RELATIVE MACHINE PRECISION.
C      THIS FUNCTION SUBPROGRAM IS A SLIGHTLY MODIFIED TRANSLATION OF
C  THE ALGOL 60 PROCEDURE  ZERO  GIVEN IN  RICHARD BRENT, ALGORITHMS
C  FOR MINIMIZATION WITHOUT DERIVATIVES, PRENTICE-HALL, INC. (1973).
C
      REAL  A,B,C,D,E,EPS,FA,FB,FC,TOL,XM,P,Q,R,S,ATOL,RTOL
      REAL SPMPAR
C
C  COMPUTE EPS, THE RELATIVE MACHINE PRECISION
C
      EPS = SPMPAR(1)
C
C INITIALIZATION
C
      A = AX
      B = BX
      FA = F(A)
      FB = F(B)
      ATOL = 0.5*AERR
      RTOL = AMAX1(0.5*RERR,2.0*EPS)
C
C BEGIN STEP
C
   10 C = A
      FC = FA
      D = B - A
      E = D
   20 IF (ABS(FC) .GE. ABS(FB)) GO TO 40
      A = B
      B = C
      C = A
      FA = FB
      FB = FC
      FC = FA
C
C CONVERGENCE TEST
C
   40 TOL = RTOL*ABS(B) + ATOL
      XM = 0.5*(C - B)
      IF (ABS(XM) .LE. TOL) GO TO 90
      IF (FB .EQ. 0.0) GO TO 90
C
C IS BISECTION NECESSARY
C
      IF (ABS(E) .LT. TOL) GO TO 70
      IF (ABS(FA) .LE. ABS(FB)) GO TO 70
C
C IS QUADRATIC INTERPOLATION POSSIBLE
C
      IF (A .NE. C) GO TO 50
C
C LINEAR INTERPOLATION
C
      S = FB/FA
      P = 2.0*XM*S
      Q = 1.0 - S
      GO TO 60
C
C INVERSE QUADRATIC INTERPOLATION
C
   50 Q = FA/FC
      R = FB/FC
      S = FB/FA
      P = S*(2.0*XM*Q*(Q - R) - (B - A)*(R - 1.0))
      Q = (Q - 1.0)*(R - 1.0)*(S - 1.0)
C
C ADJUST SIGNS
C
   60 IF (P .GT. 0.0) Q = -Q
      P = ABS(P)
C
C IS INTERPOLATION ACCEPTABLE
C
      IF ((2.0*P) .GE. (3.0*XM*Q - ABS(TOL*Q))) GO TO 70
      IF (P .GE. ABS(0.5*E*Q)) GO TO 70
      E = D
      D = P/Q
      GO TO 80
C
C BISECTION
C
   70 D = XM
      E = D
C
C COMPLETE STEP
C
   80 A = B
      FA = FB
      IF (ABS(D) .GT. TOL) B = B + D
      IF (ABS(D) .LE. TOL) B = B + SIGN(TOL,XM)
      FB = F(B)
      IF ((FB*(FC/ABS(FC))) .GT. 0.0) GO TO 10
      GO TO 20
C
C DONE
C
   90 ZEROIN = B
      RETURN
      END
