I was wondering if anyone wanted to give a shot at finding the eigenvalues of this $4\times 4$ system. I have tried without success, mainly because I end up having to solve a very nasty cubic polynomial that wolfram seems to be having trouble with as well. There may be some method to solve this that I didn't try, so I wanted to post this here just in case someone could help me out, $$J(\epsilon)=\begin{bmatrix} -\mu &&0&&0 && -\beta b/\mu\\ 0&& -\sigma && 0&& \beta b/\mu\\0&& \sigma&& -(\phi+d)&& 0\\ 0 &&0&&\alpha&&-(\phi+\epsilon)\end {bmatrix}$$
Thanks!
After removing the obvious eigenvalue $-\mu$, you are indeed left to solve the cubic polynomial
$$ {\lambda}^{3}\mu+ \left( d\mu+\epsilon\,\mu+2\,\mu\,\phi+\mu\,\sigma \right) {\lambda}^{2}+ \left( d\epsilon\,\mu+d\mu\,\phi+d\mu\,\sigma+ \epsilon\,\mu\,\phi+\epsilon\,\mu\,\sigma+\mu\,{\phi}^{2}+2\,\mu\,\phi \,\sigma \right) \lambda-\alpha\,\sigma\,\beta\,b+d\epsilon\,\mu\, \sigma+d\mu\,\phi\,\sigma+\epsilon\,\mu\,\phi\,\sigma+\mu\,{\phi}^{2} \sigma =0$$ which seems to have no special algebraic features: in fact, for suitable values of the parameters I believe it can be any cubic polynomial. A "closed-form" expression in radicals expressed in terms of all these parameters will be a nightmare. In practice, for particular numerical values of the parameters you may be best off using a numerical solver, but you can also use the general formula.