>>> import sympy
>>> import sympy.solvers as sol

>>> x = sympy.symbols('x')
>>> sol.solve(x**3-2003*x**2+6002*x-4000,x)
[1, 2, 2000]