Sympy求解方程组如何得到指定自变量表示的目标变量表达式

求解多变量方程组时，期望得到形如y = f(z)的解，但使用Sympy求解时无法得到完全消去中间变量的结果。
复现问题的简易示例代码如下：

import sympy
x = sympy.symbols("x", real=True)
y = sympy.symbols("y", real=True)
z = sympy.symbols("z", real=True)

eq0 = y - x**2
eq1 = x - z + 4

system = [eq0, eq1]
sympy.solve(system, y)[y]

上述代码运行输出为y = x**2，该结果数学层面正确但不符合需求，期望得到y = f(z)形式的结果，即y = (z - 4)**2。
核心疑问：Sympy的求解器是否支持指定输入自变量与输出因变量，实现类似如下形式的调用：

sympy.solve(equations, input=z, output=y)

关于代入替换法的补充说明
上述简易示例中可以通过将中间变量x替换为z-4得到目标结果，但实际业务场景共有36个方程，手动逐次替换的实现难度极高，因此寻求无需手动替换的原生求解方案。

Sympy原生支持该需求，不需要手动逐次替换，核心逻辑是：solve函数的第二个参数为需要被消元求解的符号集合，所有传入该参数的符号，最终都会被表示为未传入该参数的剩余符号的表达式。
之前的写法只传入了y作为待求解符号，因此求解器只会把y表示为其他所有符号（x、z）的表达式，不会主动消去中间变量x。只要把所有需要消去的中间变量、目标因变量全部传入待求解符号列表，求解器就会自动完成全量消元。
以上述示例为例，正确写法如下：

import sympy
x = sympy.symbols("x", real=True)
y = sympy.symbols("y", real=True)
z = sympy.symbols("z", real=True)

eq0 = y - x**2
eq1 = x - z + 4
system = [eq0, eq1]

# 待求解列表传入所有需要消去的中间变量x + 目标因变量y
sol = sympy.solve(system, [x, y], dict=True)[0]
print(sol[y])
# 输出结果：(z - 4)**2

如果方程数量多、中间变量多，不需要手动枚举所有待消元符号，可以通过代码自动生成待求解列表，仅保留指定的输入自变量即可：

# 收集方程组中所有符号，排除指定的输入自变量z，剩下的全部作为待消元符号
solve_targets = [s for s in sympy.sympify(system).free_symbols if s != z]
# 求解后直接取目标因变量y的表达式
y_expr = sympy.solve(system, solve_targets, dict=True)[0][y]

如果方程组存在多组解，solve会返回包含多个解字典的列表，按需筛选符合约束的解即可。

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