from sympy import *
X,Y,Z = symbols("X,Y,Z")
variables_dimless = [X,Y,Z]
Rxy,Rxz, Ry,Ryx,Ryz, Rz,Rzy,Rzz = symbols("Rxy,Rxz, Ry,Ryx,Ryz, Rz,Rzy,Rzz",positive=True)
dimensionless_DES_list = [X*(1-X),
Ry *(Y* (1 - Y) + Rzy * Z *Y ),
Rz *(Z*(1-Z)+Rxz*X*Z+Ryz*Y*Z)]
equilibria = solve( dimensionless_DES_list,(variables_dimless))
Out[2]:
[(0, 0, 0),
(0, 0, 1),
(0, 1, 0),
(0, -(Rzy + 1)/(Ryz*Rzy - 1), -(Ryz + 1)/(Ryz*Rzy - 1)),
(1, 0, 0),
(1, 0, Rxz + 1),
(1, 1, 0),
(1, -(Rxz*Rzy + Rzy + 1)/(Ryz*Rzy - 1), -(Rxz + Ryz + 1)/(Ryz*Rzy - 1))]