I am trying to obtain the steady state points of the following ODE system for a SIR model (steady state points are the solutions when the equations of the system equal zero): import sympy as sp S, I, R, b, g = sp.symbols('S I R b g') dSdt = -b * S * I dIdt = b * S * I - g * I dRdt = g * I steady_states = sp.solve([dSdt, dIdt, dRdt], [S, I, R]) steady_states When I set numerical values for b and g, it shows all solutions, but when I leave them as symbols, it only shows this solution: [(S,0,R)] the other solution should be: [(g/b,0,R)] is it possible to get both of them ? Perhaps another function in sympy? Thanks in advance. Continue reading...