First and foremost, your expression of interest involves transcendental function not defined in polynomial rings. Working in suh a ring is pointless.

Second : there is no known *general* way to solve mixed-expressions (in your case, yourexpression uses polynomials, trig and exp functions) in *closed form*. Trying to find numerical approximations of your roots is useful, but will limit your conclusions to the interval you're exploring.

Third, from `find_root?`

(emphases mine) :

Numerically find **a** root of "f" on the closed interval [a,b] (or [b,a]) if possible, where "f" is a function in the **one** variable. Note: this function only works in fixed (machine) precision, it is not possible to get arbitrary precision approximations with it.

What you neeed is a function *recursively* exploring your search interval for numerical roots and returning the (possibly empty) list of such zeros. One possible implementation :

```
def find_roots(ex, a, b):
"""
(Recursively) numerically explore the interval [a, b] fpr zeroes of
the single variable expression f, and return the (posibly empty) list
of such zeroes.
This will never return if this list is infinite
(e. g. sin(1/x) between 0 and 1).
"""
if len(ex.variables())!=1 : raise RuntimeError(
"""
find_roots :
the expression to solve must have exactly one variable
""")
if (a>b): a, b=b, a
epsilon=1e-4
try:
r1 = find_root(ex, a, b)
return find_roots(ex, a, r1-epsilon) +\
[r1] +\
find_roots(ex, r1+epsilon, b)
except RuntimeError:
return []
```

**NOTE :** The timings reported below are abnormal, and I can't reproduce them : running the same code in a newly-started Sage is *much* faster. Discard the timings and timing-related temarks...

Let's put this to work : zeroes :

```
sage: %time Zeroes = find_roots(f, -5, 5) ; Zeroes
CPU times: user 34.4 s, sys: 24 ms, total: 34.4 s
Wall time: 34.4 s
[-4.8863727269823904,
-2.248189884926029,
1.591391293749397,
2.9760183676051764]
```

Extrema :

```
sage: %time Extrema = find_roots(f.diff(x), -5, 5)
CPU times: user 26.4 s, sys: 8.01 ms, total: 26.4 s
Wall time: 26.4 s
sage: %time Maxima = [u for u in Extrema if f.diff(x,2)(x=u)<0] ; Maxima
CPU times: user 8.55 s, sys: 12 ms, total: 8.56 s
Wall time: 8.57 s
[-4.120692958580293, 2.369318406384495]
sage: %time Minima = [u for u in Extrema if f.diff(x,2)(x=u)>0] ; Minima
CPU times: user 8.43 s, sys: 6 µs, total: 8.43 s
Wall time: 8.43 s
[-1.167409895786979]
```

Check that we have no pathological cases (simultaneous nullity of f' and f'') :

```
sage: set(Extrema)-(set(Maxima)|set(Minima))
set()
```

Inflection points :

```
sage: %time Inflections = find_roots(f.diff(x, 2), -5, 5) ; Inflections
CPU times: user 1min 41s, sys: 23 µs, total: 1min 41s
Wall time: 1min 41s
[-3.340415675479655,
-0.7259000686402566,
-0.4752900420219365,
1.4375916288836839,
3.997537922882414]
```

These results appear consistent with the plots on the interval :

One notes that this search isn't especially fast. You may try to define `f`

, `f'`

and `f''`

as Cython functions : with appropriate type declaration and unique computation of the numerical constants, this should give you a notable speed boost.

Similarly, it might be interesting to reimplement `find_roots`

as a generator (with some interesting edge cases...).

All of this is left to you as an exercise (and an exercise in laziness for me :-)).

Reformatted for legibility.