I need to implement a root-searching algorithm for a function, that can possibly have some (non-removable) divergence points (DPs), e.g. f(x) = x**2 * np.tan(x).

The problem is, that the usual condition

f(x1)*f(x2) < 0 signifies a root within the interval [x1, x2]

cannot be applied on discontinuous functions. If done so on the aforementioned function, the algorithm identifies the DPs as roots.

I tried to distinguish the DPs from the roots by this criterion:

# assume x1, x2 as interval points identified by the f(x1)*f(x2) < 0 condition
x3 = (x1+x2)/2
f1, f2, f3 = f(x1), f(x2), f(x3)
if abs(f1-f2) < (abs(f1-f3)+abs(f2-f3))*0.9999999999999:
    print("DP possibly found")

which is kind of cumbersome to me (and probably also not correct, as it depends on the precision of the calculation a lot). I put there the 0.9999999999999 just to avoid possible numeric errors, but sometimes the difference between the operands can be very tiny.

(The idea of that criterion was that if the function is non-monotonous, there is probably a DP. The monotonity check is just that the distance of function values at the edges should be the same as the sum of distances of function values at left&center and center&right. I am aware of some (pathologic) examples that could mis-identify a DP, but I thought that it should work.)

So, what is the correct/better way to identify DPs?

I am looking for a solution without external modules (e.g. sympy) if possible.

I posted my testing code as a MWE here (its a bit longer mainly due to docstrings).

There is a similar question, but nobody seems to have answered well there (the accepted answer is wrong):
Finding roots but not asymptotes of a function