I am currently working on a spatially explicit ODE model with dispersion to study the population dynamics of mosquitoes. I wish to compute the equilibrium values of the populations as functions of the parameters in the model. The set of equations are presented below. $$ \begin{align} \dfrac{dE_p}{dt} =& \beta_p M_p - \dfrac{E_p}{d_{E_{p}}} - \mu^0_{E_{p}} \left( 1 + \dfrac{E_p + L_p}{K_p}\right) E_p, \\[5pt] \dfrac{dL_p}{dt} =& \dfrac{E_p}{d_{E_{p}}} - \dfrac{L_p}{d_{L_{p}}} - \mu^0_{L_{p}} \left( 1 + \gamma \dfrac{E_p + L_p}{K_p}\right) L_p, \\[5pt] \dfrac{dP_p}{dt} =& \dfrac{L_p}{d_{L_{p}}} - \dfrac{P_p}{d{_{P_{p}}}} - \mu_{P_{p}} P_p, \\[5pt] \dfrac{dM_p}{dt} =& \dfrac{1}{2} \dfrac{P_p}{d{_{P_{p}}}} - \mu_{M_{p}} M_p + \sum_{i=1}^{\overline{p}} m_{ip} M_i - \sum_{i=1}^{\overline{p}} m_{pi} M_p\\ \end{align} $$ I have equated the right hand sides to zero and tried to solve the system by susbstitution, but failed to do it. solve from Sagemath is also unable to help.
I will appreciate any help in finding a symbolic solution to this system of equations.