你遇到的情况是因为SymPy的solve()函数优先返回精确解析解，而教材里给出的是近似实数值解——这些解需要通过代数拆分+数值求解的方式才能得到，下面一步步拆解正确的求解过程：

1. 拆分方程组，分情况讨论 #

先把第一个方程因式分解：

eq1 = 20*x*y -10*x -4*x**3 = 0 → x*(20y - 10 -4x²) = 0

这自然分成两种情况：

• 情况1：x=0
  代入第二个方程eq2 =10x² -8y -8y³=0，得到：

  -8y -8y³=0 → -8y(1+y²)=0
  

  解得y=0或y=±i（复数解），这就是SymPy返回的(0,0),(0,-i),(0,i)，其中只有(0,0)是有实际意义的实数值解。

• 情况2：x≠0
  此时可以从eq1中解出y的表达式，把x作为自变量：

  20y = 10 +4x² → y = (2x² +5)/10
  

2. 代入第二个方程，转化为一元高次方程 #

把y=(2x² +5)/10代入eq2，会得到一个关于x的6次多项式方程：

10x² -8*((2x²+5)/10) -8*((2x²+5)/10)³ = 0

这个方程没有简洁的解析解，无法用solve()直接得到精确表达式，需要用数值求解工具计算实根。

3. 用SymPy的nsolve()求数值解 #

在SymPy中使用nsolve()来寻找这个方程的实根（因为方程是偶函数，负根是正根的相反数，只需要求正根即可）：

from sympy import symbols, nsolve, Eq

x, y = symbols('x y')
# 定义原方程组
eq1 = Eq(20*x*y -10*x -4*x**3, 0)
eq2 = Eq(10*x**2 -8*y -8*y**3, 0)

# 从eq1解出x≠0时的y表达式
y_expr = solve(eq1, y)[0]
# 代入eq2得到关于x的一元方程
substituted_eq = eq2.subs(y, y_expr).simplify()

# 求两个正实根（设置不同初始值）
x1 = nsolve(substituted_eq, x, 2)  # 初始值设为2，对应较大的根
y1 = y_expr.subs(x, x1)
x2 = nsolve(substituted_eq, x, 0.8)  # 初始值设为0.8，对应较小的根
y2 = y_expr.subs(x, x2)

# 输出近似结果
print(f"(±{round(float(x1),2)}, {round(float(y1),2)})")
print(f"(±{round(float(x2),2)}, {round(float(y2),2)})")

运行后就能得到教材里的结果：(±2.64, 1.90)和(±0.86, 0.65)。

• SymPy的solve()更擅长寻找有解析形式的精确解，但遇到高次多项式这类无简洁解析解的情况，必须用nsolve()这类数值求解工具。
• 对于多元方程组，先通过因式分解拆分情况，能大幅简化求解过程，避免遗漏实根。

内容的提问来源于stack exchange，提问作者Munseok