The answer from Maelstrom is good but I just want to add a few points.

The values you substitute are all floats and with those values the polynomial is ill-conditioned. That means that the form of the expression that you substitute into can affect the accuracy of the returned results. That is one reason why substituting values into the solution from solve does not necessarily give exactly the same value that you get from substituting before calling solve.

Also before you substitute the symbols it isn’t possible for solve to know which of the three roots is real. That’s why you get three solutions from `solve(expr, r)` and only one solution from `solve(expr.subs(vdict), r)`. The third solution which is real after the substitution is the same (ignoring the tiny imaginary part) as returned by solve after the substitution:

``````In [7]: soln1[2].subs(vdict).n()
Out[7]: 0.000182122477993494 + 1.23259516440783e-32⋅ⅈ

In [8]: solve(expr.subs(vdict), r)
Out[8]: [0.000182122477993494]
``````

Because the polynomial is ill-conditioned and has a large gradient at the root `nsolve` has a hard time finding this root. However `nsolve` can find the root if given a narrow enough interval:

``````In [9]: nsolve(expr.subs(vdict), r, [0.0001821, 0.0001823])
Out[9]: 0.000182122477993494
``````

Since this is essentially a polynomial your best bet is actually to convert it to a polynomial and use `nroots`. The quickest way to do this is using `as_numer_denom` although in this case that introduces a spurious root at zero:

``````In [26]: Poly(expr.subs(vdict).as_numer_denom()[0], r).nroots()
Out[26]: [0, 0.000182122477993494, -9.17942953565356e-5 - 0.000158143657514284⋅ⅈ, -9.17942953565356e-5 + 0.000158143657514284⋅ⅈ]
``````