Function ZBRENT(X1, X2, TOL)
ITMAX = 100
EPS = 0.00000003
A = X1
B = X2
FA = FUNC(A)
FB = FUNC(B)
If FB * FA > 0# Then
Print "Root must be bracketed for ZBRENT"
End If
FC = FB
For ITER = 1 To ITMAX
If FB * FC > 0 Then
C = A
FC = FA
D = B - A
E = D
End If
If Abs(FC) < Abs(FB) Then
A = B
B = C
C = A
FA = FB
FB = FC
FC = FA
End If
TOL1 = 2# * EPS * Abs(B) + 0.5 * TOL
XM = 0.5 * (C - B)
If Abs(XM) ZBRENT = B
Exit Function
End If
If Abs(E) >= TOL1 And Abs(FA) > Abs(FB) Then
S = FB / FA
If A = C Then
P = 2# * XM * S
Q = 1 - S
Else
Q = FA / FC
R = FB / FC
P = S * (2# * XM * Q * (Q - R) - (B - A) * (R - 1#))
Q = (Q - 1#) * (R - 1#) * (S - 1#)
End If
If P > 0# Then Q = -Q
P = Abs(P)
If 3# * XM * Q - Abs(TOL1 * Q) < Abs(E * Q) Then
AAA = 3# * XM * Q - Abs(TOL1 * Q)
Else
AAA = Abs(E * Q)
End If
If 2# * P < AAA Then
E = D
D = P / Q
Else
D = XM
E = D
End If
Else
D = XM
E = D
End If
A = B
FA = FB
If Abs(D) > TOL1 Then
B = B + D
Else
B = B + Abs(TOL1) * Sgn(XM)
End If
FB = FUNC(B)
Next ITER
Print "ZBRENT exceeding maximum iterations."
ZBRENT = B
End Function
