1 REAL FUNCTION ZEROIN (AX,BX,F,TOL)

2 REAL AX,BX,F,TOL

3 C

4 C

5 REAL A, B, C, D, E, EPS, FA, FB, FC, TOL1, XM, P, Q, R, S

6 C

7 C COMPUTER EPS, THE RELATIVE MACHINE PRECISION

8 C

9 EPS = 1.0

10. DO

EPS = EPS/2.0

11 TOL1 = 1.0 + EPS

12 WHILE (TOL1 .GT. 1.0)

13 C

14 C INITIALIZATION

15 C

18 A = AX

19 B = BX

20 FA = F(A)

21 FB = F(B)

22 C

23 C BEGIN STEP

24 C

25 20 C = A

26 FC = FA

27 D = B - A

28 E = D

31 IF (ABS(FC) .GE. ABS(FB)) GO TO 40

32 A = B

33 B = C

34 C = A

35 FA = FB

36 FB = FC

37 FC = FA

38 C

39 C CONVERGENCE TEST

40 C

41 40 TOL1 = 2.0*EPS*ABS(B) + 0.5*TOL

42 XM = .5*(C-B)

43 IF ((ABS(XM) .LE. TOL1) .OR.

44 (FB .EQ. 0.0)) GO TO 90

45 C

46 C IS BISECTION NECESSARY

47 C

48 IF ((ABS(E) .LT. TOL1) .OR.

49 (ABS(FA) .LE. ABS(FB))) GO TO 70

50 C

51 C IS QUADRATIC INTERPOLATION POSSIBLE

52 C

53 IF (A .EQ. C)

54 C

55 C LINEAR INTERPOLATION

56 C

57 S = FB/FA

58 P = 2.0*XM*S

59 Q = 1.0 - S

60 ELSE

61 C

62 C INVERSE QUADRATIC INTERPOLATION

63 C

64 Q = FA/FC

65 R = FB/FC

66 S = FB/FA

67 P = S*(2.0*XM*Q*(Q-R) - (B-A) * (R-1.0))

68. Q = (Q-1.0)*(R-1.0)*(S-1.0)

END IF

69 C

70 C ADJUST SIGNS

71 C

72 IF (P .GT. 0.0) Q = -Q

73 P = ABS(P)

74 C

75 C IS INTERPOLATION ACCEPTABLE

76 C

77 IF (((2.0*P) .GE. (3.0*XM*Q - ABS(TOL1*Q))) .OR.

78 (P .GE. ABS(0.5*E*Q))) GO TO 70

79 E = D

80 D = P/Q

81 GO TO 80

82 C

83 C BISECTION

84 C

85 70 D = XM

86 E = D

87 C

88 C COMPLETE STEP

89 C

90 80 A = B

91 FA = FB

92. IF (ABS(D) .GT. TOL1)
93. B = B + D
94. ELSE
95. B = B + SIGN(TOL1,XM)

END IF

94 FB = F(B)

95 IF ((FB*(FC/ABS (FC))) .GT. 0.0) GO TO 20

96 GO TO 30

97 C

98 C DONE

99 C

100 90 ZEROIN = B

101 RETURN

102 END