Resolução de sistema de equações não lineares com python

Posso resolver um sistema de equações não lineares em termos de parâmetros em python? Existe um exemplo ou tutorial? Eu posso fazer isso facilmente no bordo, mas as expressões para o meu sistema particular são muito grandes e copiá-las é muito difícil.

Exemplo:

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

Você deve obter as soluções:

[[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]]

A razão pela qual pergunto: tenho um grande sistema de ODEs não lineares. Quero resolver os pontos fixos (isso é possível, já foi feito no bordo, mas eles são grandes e feios). Quero criar mais expressões a partir dos pontos fixos e depois usar o pacote de otimização no scipy. Prefiro fazer tudo em python do que traduzir as coisas de um lado para o outro, uma vez que é muito ineficiente e podem ser cometidos erros.

questionAnswers(3)

yourAnswerToTheQuestion