Sistema de resolución de ecuaciones no lineales con python

¿Puedo resolver un sistema de ecuaciones no lineales en términos de parámetros en python? ¿Hay un ejemplo o tutorial? Puedo hacer esto fácilmente en arce, pero las expresiones para mi sistema en particular son bastante grandes y copiarlas es bastante difícil.

Ejemplo:

sigma*(y-x) = 0
x*(rho-z)-y = 0
x*y-beta*z = 0

Deberías obtener las soluciones:

[[x = 0, y = 0, z = 0], [x = sqrt(beta*rho-beta), y = sqrt(beta*rho-beta), z = rho-1],
[x = -sqrt(beta*rho-beta), y = -sqrt(beta*rho-beta), z = rho-1]]

La razón por la que pregunto: tengo un gran sistema de EDO no lineales. Quiero resolver los puntos fijos (esto es factible, se ha hecho en arce, pero son grandes y feos). Quiero crear más expresiones a partir de los puntos fijos y luego usar el paquete de optimización en scipy. Prefiero hacerlo todo en Python que traducir cosas de un lado a otro, ya que es muy ineficiente y se pueden cometer errores.

Respuestas a la pregunta(3)

Su respuesta a la pregunta