How to get all the zeroes for the derivative of a function using DualNumbers

Question

The original problem I had was, given an interval, [a, b] and a function f(x) and its derivative continuous over that interval, find the minimum and maximum of the derivative over that integral, (and since it's continuous, every single value of the derivative will be within that interval), and I'd have to do this using DualNumber.

So I implemented DualNumbers in python and I was able to operator overload in order to achieve the functioning of the class, and I also implemented the function which merely calculated the derivative of the function at a specific point.

I wasn't sure how would I find the minimum and maximum of the derivative since it would include real numbers as well, and hence I set out to find the zero's of the derivative (do note there might be multiple 0's in the interval).

I wasn't sure how I would go about that, and hence I'm here asking this question.

Here is my DualNumbers class implementation:

class DualNumbers:
    def __init__(self, coef = 0, epsilon = 0) -> None:
        self.coef = coef
        self.epsilon = epsilon

    def __add__(self, __o):
        if type(__o) != DualNumbers:
            if type(__o) == float or type(__o) == int:
                return DualNumbers(self.coef + __o, self.epsilon)
            else:
                raise "Incorrect variable type called in add"
        
        return DualNumbers(self.coef + __o.coef, self.epsilon + __o.epsilon)

    def __radd__(self, __o):
        if type(__o) != DualNumbers:
            if type(__o) == float or type(__o) == int:
                return DualNumbers(self.coef + __o, self.epsilon)
            else:
                raise "Incorrect variable type called in add"
        
        return DualNumbers(self.coef + __o.coef, self.epsilon + __o.epsilon)

    def __sub__(self, __o):
        if type(__o) != DualNumbers:
            if type(__o) == float or type(__o) == int:
                return DualNumbers(self.coef - __o, self.epsilon)
            else:
                raise "Incorrect variable type called in add"
        
        return DualNumbers(self.coef - __o.coef, self.epsilon - __o.epsilon)

    def __rsub__(self, __o):
        if type(__o) != DualNumbers:
            if type(__o) == float or type(__o) == int:
                return DualNumbers(__o - self.coef, -1 * self.epsilon)
            else:
                raise "Incorrect variable type called in add"
        
        return DualNumbers(__o - self.coef, -1 * self.epsilon)

    def __eq__(self, __o: object) -> bool:
        if type(__o) != DualNumbers:
            return False
        
        return self.coef == __o.coef and self.epsilon == __o.epsilon

    def __ne__(self, __o: object) -> bool:
        return not self.__eq__(__o)

    def __mul__(self, __o) -> float:
        if type(__o) != DualNumbers:
            if type(__o) == int or type(__o) == float:
                return DualNumbers(self.coef * __o, self.epsilon * __o)
            else:
                raise TypeError
        
        return DualNumbers(self.coef * __o.coef, (self.epsilon * __o.coef) + (self.coef * __o.epsilon))
    
    def __rmul__(self, __o) -> float:
        if type(__o) != DualNumbers:
            if type(__o) == int or type(__o) == float:
                return DualNumbers(self.coef * __o, self.epsilon * __o)
            else:
                raise TypeError
        
        return DualNumbers(self.coef * __o.coef, (self.epsilon * __o.coef) + (self.coef * __o.epsilon))

    def __div__(self, __o: object) -> float:
        if type(__o) != DualNumbers:
            return False
        
        try:
            return DualNumbers(self.coef / __o.coef, ((self.epsilon * __o.coef) - (__o.epsilon * self.coef)) / (__o.coef ** 2))
        except:
            raise "c is 0, division with 0"

    def __str__(self) -> str:
        return (f"{self.coef} + {self.epsilon}ε")

    def __repr__(self) -> str:
        return (f"{self.coef} + {self.epsilon}ε")

    def __pow__(self, other):
        x = DualNumbers(self.coef, self.epsilon)
        a = x
        for _ in range(other - 1):
            x = x * a
        return x

and my derivative computing function is:

def compute_derivative(f, x):
    x_dn = DualNumbers(x, 1)
    
    y_dn = f(x_dn)
    
    derivative = y_dn.epsilon
    
    return derivative

Is there a better way to do this, or worst case could I use the second derivative (not sure how I'd go about calculating that) to get the 0's in some way (Newton's method or something)?