SymPy is giving me complex solutions, even when I constrain my variables to be real, as in this answer.
Here's an example finding the intersection of a unit sphere and the plane X=1. There should be exactly 1 real solution, (1,0,0).
>>> import sympy
>>> x,y,z = sympy.symbols('x,y,z', real=True)
>>> list(sympy.solve([x**2 + y**2 + z**2 - 1, x - 1], [x,y,z]))
[(1, -I*z, z), (1, I*z, z)]
These are correct, but not real.
I can't use solveset(…, domain=S.Reals)
because solveset doesn't seem to support systems of equations.