先指出一个错误：你list1中只有一个元素，应该用list1[0]取出；
为留坝等地区用户提供了全套网页设计制作服务，及留坝网站建设行业解决方案。主营业务为网站制作、网站建设、留坝网站设计，以传统方式定制建设网站，并提供域名空间备案等一条龙服务，秉承以专业、用心的态度为用户提供真诚的服务。我们深信只要达到每一位用户的要求，就会得到认可，从而选择与我们长期合作。这样，我们也可以走得更远！
这个问题，我也一直在找解决办法，苦搜无果，自己想到了增加方程组变量的方法来解决：新增变量-表达式=0，把这个方程同之前你得到的结果组成三元一次方程组，得出新增变量的解即可。
z=Symbol('z')
result1=solve([z-list1[0],x-result[x],y-result[y]],[x,y,z])
result1[z]就是你要的结果，拿走不谢~~
如果一元方程的解（x）带回表达式，思路一样，只是注意一元方程的解是存放在列表里（假设为result[]），而不是字典，列表中的第一个元素为实数解，所以代码变为：
y=Symbol('y')
result1=solve([y-list1[0],x-result[0],[x,y])
result1[y]即是。
用牛顿迭代法 + 多项式除法化简。
1）针对方程组 f(x)，首先用牛顿迭代法得到方程的第一个根（a），那么 f（x） = (x-a)g(x)
2）用多项式除法，计算 g(x) = f(x)/(x-a)
重复第一步，得到 g(x) 的根，然后再重复第二步，进一步对方程降幂。
最终就可以化简整个方程。
不是很明确你需要做到什么程度，但基本可以通过以下两个手段得到：
手工解方程得到解析解，然后套入公式
使用一些工具包例如numpy可以自动求解
以下都给出例子
import numpy as np
import matplotlib.pyplot as plt
plt.axis("equal")
a = np.linspace(1,10,100) # a 的变化范围可以自己挑，前两个参数控制，
# 使用 numpy 自动求解
res = []
for x in a:
A = np.mat("1, 2; {}, -1".format(x))
b = np.mat("{}, 10".format(x)).T
res.append(np.linalg.solve(A, b))
# 计算完毕后取出每对x和y
x1 = [float(r[0]) for r in res]
y1 = [float(r[1]) for r in res]
plt.plot(x1, y1)
#####################################
# 手工计算过程很简单不放上来了，直接上结果
x2 = [(a1 + 20) / (2*a1 + 1) for a1 in a]
y2 = [(a1**2 - 10) / (2*a1 + 1) for a1 in a]
plt.plot(x2, y2)