I'm messing around with writing some SymPy code to handle symbolic expressions with imaginary numbers.

To start out, I want to get it to take x and y as real numbers and find the solution where x=iy. So I can do this as follows.

x, y = sympy.symbols("x y", real=True)  
print(sympy.solve([x-sympy.I*y]))

(SymPy solve takes a list of values, all of which must be 0. So x-iy=0 => x=iy). SymPy will correctly tell me

[{x: 0, y: 0}]

However, if I do this a (theoretically identical) way:

x, y = sympy.symbols("x y")
print(sympy.solve([x-sympy.I*y, sympy.im(y), sympy.im(x)]))

Then now SymPy tells me

[{re(y): y, re(x): I*y, im(x): 0, x: I*y, im(y): 0}]

And this is technically correct, but hasn't done everything for me. Is this just a limitation in SymPy, or can I get it to give me x=y=0 by constraining complex x and y in this way?

Because SymPy is better at simplifying pairs of real numbers than complex numbers, the following strategy helps: set up real variables for real/imaginary parts, then form complex variables from them.

from sympy import *
x1, x2, y1, y2 = symbols("x1 x2 y1 y2", real=True)  
x = x1 + I*x2
y = y1 + I*y2

Now x and y can be used as complex variables in an equation such as yours

sol = solve([x-I*y, im(y), im(x)])
print(x.subs(sol[0]), y.subs(sol[0])) 

Output: 0 0.