# Routh-Hurwitz criterion not giving correct answer when done manually?

Consider the system:

\begin{align} \frac{dS}{dt} &= \nu N -\frac{\beta S I}{N} + \xi R - \nu S\\ \frac{dE}{dt} &= \frac{\beta S I}{N}- \sigma E -\nu E \\[2ex] \frac{dI}{dt} &= \sigma E -\gamma I -\nu I \\[2ex] \frac{dR}{dt} &= \gamma I -\xi R - \nu R \end{align}

The Jacobian is: \begin{align} J\left(S,E,I\right) = \begin{bmatrix} -\frac{\beta I_2^*}{N}-B &-\xi & -\frac{\beta S_2^*}{N}-\xi \\[1ex] \frac{\beta I_2^*}{N} & -C & \frac{\beta S_2^*}{N}\\[1ex] 0 & \sigma & -D \end{bmatrix}. \end{align} where $$S_2^*=\frac{CDN}{\beta \sigma}$$ $$I_2^* = \frac{BN(-C D + \beta \sigma)}{\beta (C D + D \xi + \xi \sigma)}$$

Where $$B = \xi + \nu, \quad C = \sigma+\nu \quad \text{and } D = \gamma + \nu$$.

\begin{align*} P\left(\lambda\right) &= \lambda ^3 + \left[B +C +D + \frac{B\left(\beta \sigma -CD \right)}{CD + \xi\left(D+\sigma\right)} \right]\lambda^2\\[2ex] & \quad+\left[BC + BD + \frac{BC\left(\beta \sigma -CD \right)}{CD + \xi\left(D+\sigma\right)} +\frac{BD\left(\beta \sigma -CD \right)}{CD + \xi\left(D+\sigma\right)}+\frac{B\xi \left(\beta \sigma -CD \right)}{CD + \xi\left(D+\sigma\right)} \right]\lambda\\[2ex] & \quad+\left[B\left(\beta \sigma -CD\right) \right] \end{align*}

Using Chris K's EcoEvo package to test the stability we see:

Simplify[EcoStableQ[eq2[[2]]]]


outputs:

True when $$CD<\beta \sigma$$ and false when $$CD>\beta \sigma$$, which is as expected.

Now when I try to do it with brute force using the Routh-Hurwitz conditions:

\begin{align} a_1 & > 0\\[1ex] a_3 & > 0\\[1ex] a_1 a_2 & > a_3 \end{align} For our system: \begin{align*} a_1 & = B +C +D + \frac{B\left(\beta \sigma -CD \right)}{CD + \xi\left(D+\sigma\right)} \\[2.5ex]& = \frac{C^2 D + C D^2 + B D \xi + C D \xi + D^2 \xi + B \beta \sigma + B \xi \sigma + C \xi \sigma + D \xi \sigma}{CD + \xi\left(D+\sigma\right)} \\[3ex] a_2 & = BC + BD + \frac{BC\left(\beta \sigma -CD \right)}{CD + \xi\left(D+\sigma\right)} +\frac{BD\left(\beta \sigma -CD \right)}{CD + \xi\left(D+\sigma\right)}+\frac{B\xi \left(\beta \sigma -CD \right)}{CD + \xi\left(D+\sigma\right)} \\[2.5ex] & = \frac{B D^2 \xi + B C \beta \sigma + B D \beta \sigma + B C \xi \sigma + B D \xi \sigma + B \beta \xi \sigma}{CD + \xi\left(D+\sigma\right)}\\[3ex] a_3 & =B\left(\beta \sigma -CD\right) \end{align*}

Clearly $$a_1>0$$ and $$a_3>0$$ if $$\frac{\beta \sigma}{CD} >1$$ but why is the third condition failing when it should hold as we see from Chris's package?

We know a priori if $$\frac{\beta \sigma}{CD} >1$$ then the eigenvalue(s) must have negative real parts implying the equilibrium point is stable.

EDIT:

Here is the system is Mathematica(you will need to load the package in Mathematica first; PacletInstall["EcoEvo", "Site" -> "http://raw.githubusercontent.com/cklausme/EcoEvo/master"] and << EcoEvo:

r := n - s - e - i;
SetModel[{Pop[
pop] -> {Component[
s] -> {Equation :>
B n - ξ e - ξ i - (β/n) s i  - B s},
Component[e] -> {Equation :>  (β/n) s i - C e },
Component[i] -> {Equation :> σ e - D i }},
Parameters :> {β > 0, A > 0, ξ > 0, σ > 0, B > 0,
C > 0, D > 0}}]
eq2 = SolveEcoEq[]


EDIT 2:

The paper:

I believe the author has made a typo in equation (75) on the second matrix.

Their model is slightly different:

r := 1 - s - e - i;
SetModel[{Pop[
pop] -> {Component[
s] -> {Equation :>
A - \[Alpha] e - \[Alpha] i - \[Beta] s i  - B s},
Component[e] -> {Equation :>  \[Beta] s i - C e },
Component[i] -> {Equation :> \[Sigma] e - D i }},
Parameters :> {\[Beta] > 0, A > 0, \[Xi] > 0, \[Sigma] > 0, B > 0,
C > 0, D > 0}}]
eq3 = SolveEcoEq[]


Their $$S_2^*=\frac{CD}{\beta \sigma}$$, $$I_2^* = \frac{-B C D + \beta \sigma A}{\beta (C D + D \alpha + \alpha \sigma)}$$

By remarkable chance, they managed to show R-H criteria working as expected with their typo!

• What happened to the fourth variable R? Should that not also be accounted for in the Jacobian? Commented Sep 28, 2021 at 13:53
• We reduced the system since we have the equality $S+E+I+R=N$
– Math
Commented Sep 28, 2021 at 13:56

I can't say why your approach didn't work, but my RouthHurwitzCriteria function uses a simplified test for 3x3 matrices due to Fuller (1968), which I first learned about from Gandolfo (1997):

There are extensions to higher-order systems in Murata (1977) that I thought about implementing, but never had the time. Any matrix lovers out there want to help?

Fuller AT. 1968. Conditions for a matrix to have only characteristic roots with negative real parts. Journal of Mathematical Analysis and Applications 23: 71-98.

• Thank you for the references. Since you have answered a few of my questions and hence, are well aware of these models, shouldn't the R-H conditions presented in my question hold $\textit{if}$ $\mathcal{R}_0>1$?
– Math
Commented Sep 28, 2021 at 13:55
• Higher-order you mean 4x4 matrix, 4D dynamic system? But that is not needed, no? Commented Jul 21, 2022 at 12:51
• @ВалерийЗаподовников Needed by OP? No. Needed by me? Not too much, but I wouldn't mind. Commented Jul 21, 2022 at 16:18
• Are there real applications for 4x4 Hurwitz stable matrices? That is not 4 degree polinomials Hurwitz stability, Hurwitz matrix always has some zero elements... Commented Jul 21, 2022 at 21:48
• @ВалерийЗаподовников It seems potentially useful for analytical results in the narrow range of cases where the resulting algebraic conditions are not too complicated. Otherwise you're probably stuck with numerical results, in which case calculating the eigenvalues isn't such a big deal. All of which is why I haven't bothered to implement it! Commented Jul 21, 2022 at 22:06

Your approach didn't work because the last condition cannot be solved easily, that is, you need at least one parameter of degree one to solve without complexity. Then, Chris's routine does work perfectly.

However, it is possible to analyze stability if we think about the Hopf bifurcation, and we solve the following equation

$$\label{1} \left| \begin{array}{ccc} J_{11}+J_{22} & J_{23} & -J_{13} \\ J_{32} & J_{11}+J_{33} & J_{12} \\ -J_{31} & J_{21} & J_{22}+J_{33} \\ \end{array} \right|=0,\tag 1$$ for some parameter of your system. Firstly, if we consider $$\beta=\displaystyle\frac{cd}{\sigma}$$, we obtain the following three real eigenvalues:

$$\lambda_{1}=-b,\hspace{0.1cm}\lambda_{2}=-c,\hspace{0.1cm}\lambda_{3}=-d.$$ Secondly, as an example, for $$\beta>\displaystyle\frac{cd}{\sigma}$$ we set $$d=\xi=1$$, $$c=1/100$$, $$\sigma=13$$, $$\beta=\displaystyle\frac{101}{1300}$$ and solve $$(\ref{1})$$ for $$b$$.

S2 = (c d n)/(β σ); I2 = (b n (β σ - c d))/(β (c d + d ξ + ξ σ));
J = {{-((β I2)/n) - b, -ξ, -((β S2)/n) - ξ}, {(β I2)/n, -c, (β I2)/n}, {0, σ, -d}};
M = {{J[[2, 2]] + J[[1, 1]], J[[2, 3]], -J[[1, 3]]}, {J[[3, 2]],
J[[3, 3]] + J[[1, 1]], J[[1, 2]]}, {-J[[3, 1]], J[[2, 1]], J[[3, 3]] + J[[2, 2]]}};


The Hopf bifurcation values:

detM[{b_, c_, d_, β_, σ_, ξ_}] := Evaluate[Det[M]]
b1=Solve[Factor[detM[{b, 1/100, 1, 101/1300, 13, 1}]] == 0, b][[1, 1, 2]]
b2=Solve[Factor[detM[{b, 1/100, 1, 101/1300, 13, 1}]] == 0, b][[2, 1, 2]]


The transversality condition for $$b_{1}$$ and $$b_{2}$$:

Factor[-D[detM[{b, 1/100, 1, 101/1300, 13, 1}], b] /. b -> b1]
(*-(Sqrt[107236043973001]/14010000)*)
(*locally stable for b<b1*)
Factor[-D[detM[{b, 1/100, 1, 101/1300, 13, 1}], b] /. b -> b2]
(*(Sqrt[107236043973001]/14010000)*)
(*locally stable for b>b2*)


You can verify with the roots.

Appendix

Typo:

 J = {{-((β I2)/n) - b, -ξ, -((β S2)/n) - ξ}, {(β I2)/n, -c, (β I2)/n}, {0, σ, -d}};(*incorrect*)
J = {{-((β I2)/n) - b, -ξ, -((β S2)/n) - ξ}, {(β I2)/n, -c, (β S2)/n}, {0, σ, -d}}; (*correct*)


Example: If we set $$\beta=\displaystyle\frac{\beta_{0}+cd}{\sigma}$$ ($$\beta_{0}>0$$), $$d=\xi$$ and $$b=\sigma$$, your equilibrium point is locally stable.

pol = a0 λ^3 + a1 λ^2 + a2 λ + a3;


where

a0 = 1;
a1 = c + ξ + σ + (β0 σ)/(ξ (c + ξ + σ));
a2 = (σ (c^2 ξ + ξ (2 β0 + ξ (ξ + σ)) +c (β0 + ξ (2 ξ + σ))))/(ξ (c + ξ + σ)):
a3 = -β0 σ + (σ (c + ξ + σ + (β0 σ)/(ξ (c + ξ + σ))) (c^2 ξ + ξ
(2 β0 + ξ (ξ + σ)) + c (β0 + ξ (2 ξ + σ))))/(ξ (c + ξ + σ))
H2 = a1 a2 - a3;


Checking stability:

Reduce[a1 > 0 && a2 > 0 && H2 > 0 && β0 > 0 && σ > 0 && ξ > 0 && c > 0]
(*σ > 0 && ξ > 0 && β0 > 0 && c > 0*)


Stability via Hopf bifurcation:

Example: If we set $$\beta=\displaystyle\frac{\beta_{0}+cd}{\sigma}$$ we can verify the stability as follows:

 pol = a0 λ^3 + a1 λ^2 + a2 λ + a3; (*characteristic polynomial*)


where

a0 = 1;
a1 = b + c + d + (b β0)/(c d + ξ (d + σ));
a2 = (b (c^2 d + d^2 ξ + β0 ξ + d (β0 + ξ σ) + c (d^2 + β0 + d ξ
+ ξ σ)))/(c d + ξ (d + σ));
a3 = b β0;
H2 = a1 a2 - a3;
Collect[H2, b, FullSimplify]
(*(b^2 (c d + β0 + ξ (d + σ)) (c^2 d + d^2 ξ + β0 ξ + d (β0 + ξ σ) + c (β0 + d (d + ξ) + ξ σ)))/(c d + ξ (d + σ))^2 + b (-β0 + (c + d) (c + d + (β0 (c + d + ξ))/(c d + ξ (d + σ))))*)


Solve $$H_{2}=0$$ for $$b$$ and find the critical value of the Hopf bifurcation $$b=b_{0}$$. Then, analyze the stability with respect to $$b_{0}$$.

Factor[Reduce[a1 > 0 && a2 > 0 && H2 > 0 && c > 0 && d > 0 && ξ > 0 && σ > 0 && β0 > 0 && b > 0, b]]
(*Conditions "arbitrary values"*)


This is the strategy I followed. Check it out.

Remark

Symbolically, the Routh-Hurwitz test consists of solving inequalities that can be simple or very complicated. There is no loss of generality if you choose to solve the inequalities for some parameter. We only determine the conditions (inequalities) on the parameters to obtain the local stability of our equilibrium.

Remark

The following matrix:

$$A=\left( \begin{array}{ccc} a_{11}+a_{22} & a_{23} & -a_{13} \\ a_{32} & a_{11}+a_{33} & a_{12} \\ -a_{31} & a_{21} & a_{22}+a_{33} \\ \end{array} \right),$$ corresponds to the bialternate matrix product $$2 A \odot I_{3}$$. You can use a bialternate matrix product to extend the stability criteria for higher dimensions of the phase space. For more details, see Elements of Applied Bifurcation Theory - Yuri A. Kuznetsov (Appendix B).

You can use BialternateSum.

• How did you get those eigenvalues?
– Math
Commented Sep 29, 2021 at 12:38
• Can you expand on "Solve $H_2=0$ for $b$ and find the critical value of the Hopf bifurcation $b=b_0$. Then, analyze the stability with respect to $b_0$." because i don't seem to get the any answers. Strange.
– Math
Commented Oct 6, 2021 at 10:21
• Shouldn't you have solved matrix $M$ and not $J$?
– Math
Commented Oct 7, 2021 at 11:45
• @E. Chan-López: Why didn’t you check all 3 conditions of Hopf bifurcation theorem(non-hyperbolicity, transversality, genericity)?
– Leo
Commented Oct 7, 2021 at 12:16
• Can you add how to use your new FunctionRepository function? Commented Jul 21, 2022 at 12:52

Conditions 1 and 2 are trivially satisfied. Now we deal with the 3rd condition. We simplify here by multiplication of both sides with denominator and take the RHS to the LHS yielding:

B C^3 D^3 + 3 B C^2 D^3 \[Xi] + B C D^4 \[Xi] + B^2 D^3 \[Xi]^2 +
2 B C D^3 \[Xi]^2 + B D^4 \[Xi]^2 + B C^3 D \[Beta] \[Sigma] +
B C^2 D^2 \[Beta] \[Sigma] + B C D^3 \[Beta] \[Sigma] +
B C^3 D \[Xi] \[Sigma] + 4 B C^2 D^2 \[Xi] \[Sigma] +
B C D^3 \[Xi] \[Sigma] + B^2 C D \[Beta] \[Xi] \[Sigma] +
2 B C^2 D \[Beta] \[Xi] \[Sigma] + 2 B^2 D^2 \[Beta] \[Xi] \[Sigma] +
B C D^2 \[Beta] \[Xi] \[Sigma] + B D^3 \[Beta] \[Xi] \[Sigma] +
B^2 C D \[Xi]^2 \[Sigma] + B C^2 D \[Xi]^2 \[Sigma] +
2 B^2 D^2 \[Xi]^2 \[Sigma] + 5 B C D^2 \[Xi]^2 \[Sigma] +
2 B D^3 \[Xi]^2 \[Sigma] + B^2 D \[Beta] \[Xi]^2 \[Sigma] +
B C D \[Beta] \[Xi]^2 \[Sigma] + B^2 C \[Beta]^2 \[Sigma]^2 +
B^2 D \[Beta]^2 \[Sigma]^2 + 2 B^2 C \[Beta] \[Xi] \[Sigma]^2 +
B C^2 \[Beta] \[Xi] \[Sigma]^2 + 2 B^2 D \[Beta] \[Xi] \[Sigma]^2 +
B D^2 \[Beta] \[Xi] \[Sigma]^2 + B^2 \[Beta]^2 \[Xi] \[Sigma]^2 +
B^2 C \[Xi]^2 \[Sigma]^2 + B C^2 \[Xi]^2 \[Sigma]^2 +
B^2 D \[Xi]^2 \[Sigma]^2 + 3 B C D \[Xi]^2 \[Sigma]^2 +
B D^2 \[Xi]^2 \[Sigma]^2 + B^2 \[Beta] \[Xi]^2 \[Sigma]^2 +
B C \[Beta] \[Xi]^2 \[Sigma]^2 - B D \[Beta] \[Xi]^2 \[Sigma]^2 -
B \[Beta] \[Xi]^2 \[Sigma]^3 >0


Since we have $$B= \xi +\nu$$, $$C=\sigma+\nu$$ and $$D=\gamma+\nu$$, we input this into the above expression, we obtain:

\[Gamma]^3 \[Nu]^4 + 3 \[Gamma]^2 \[Nu]^5 +
3 \[Gamma] \[Nu]^6 + \[Nu]^7 + \[Gamma]^4 \[Nu]^2 \[Xi] +
8 \[Gamma]^3 \[Nu]^3 \[Xi] + 18 \[Gamma]^2 \[Nu]^4 \[Xi] +
16 \[Gamma] \[Nu]^5 \[Xi] + 5 \[Nu]^6 \[Xi] +
2 \[Gamma]^4 \[Nu] \[Xi]^2 + 14 \[Gamma]^3 \[Nu]^2 \[Xi]^2 +
30 \[Gamma]^2 \[Nu]^3 \[Xi]^2 + 26 \[Gamma] \[Nu]^4 \[Xi]^2 +
8 \[Nu]^5 \[Xi]^2 + \[Gamma]^4 \[Xi]^3 +
8 \[Gamma]^3 \[Nu] \[Xi]^3 + 18 \[Gamma]^2 \[Nu]^2 \[Xi]^3 +
16 \[Gamma] \[Nu]^3 \[Xi]^3 +
5 \[Nu]^4 \[Xi]^3 + \[Gamma]^3 \[Xi]^4 +
3 \[Gamma]^2 \[Nu] \[Xi]^4 +
3 \[Gamma] \[Nu]^2 \[Xi]^4 + \[Nu]^3 \[Xi]^4 + \[Beta] \[Gamma]^3 \
\[Nu]^2 \[Sigma] + 4 \[Beta] \[Gamma]^2 \[Nu]^3 \[Sigma] +
3 \[Gamma]^3 \[Nu]^3 \[Sigma] +
6 \[Beta] \[Gamma] \[Nu]^4 \[Sigma] +
9 \[Gamma]^2 \[Nu]^4 \[Sigma] + 3 \[Beta] \[Nu]^5 \[Sigma] +
9 \[Gamma] \[Nu]^5 \[Sigma] + 3 \[Nu]^6 \[Sigma] +
2 \[Beta] \[Gamma]^3 \[Nu] \[Xi] \[Sigma] + \[Gamma]^4 \[Nu] \[Xi] \
\[Sigma] + 10 \[Beta] \[Gamma]^2 \[Nu]^2 \[Xi] \[Sigma] +
14 \[Gamma]^3 \[Nu]^2 \[Xi] \[Sigma] +
18 \[Beta] \[Gamma] \[Nu]^3 \[Xi] \[Sigma] +
40 \[Gamma]^2 \[Nu]^3 \[Xi] \[Sigma] +
10 \[Beta] \[Nu]^4 \[Xi] \[Sigma] +
43 \[Gamma] \[Nu]^4 \[Xi] \[Sigma] +
16 \[Nu]^5 \[Xi] \[Sigma] + \[Beta] \[Gamma]^3 \[Xi]^2 \[Sigma] + \
\[Gamma]^4 \[Xi]^2 \[Sigma] +
8 \[Beta] \[Gamma]^2 \[Nu] \[Xi]^2 \[Sigma] +
15 \[Gamma]^3 \[Nu] \[Xi]^2 \[Sigma] +
19 \[Beta] \[Gamma] \[Nu]^2 \[Xi]^2 \[Sigma] +
50 \[Gamma]^2 \[Nu]^2 \[Xi]^2 \[Sigma] +
12 \[Beta] \[Nu]^3 \[Xi]^2 \[Sigma] +
62 \[Gamma] \[Nu]^3 \[Xi]^2 \[Sigma] + 26 \[Nu]^4 \[Xi]^2 \[Sigma] +
2 \[Beta] \[Gamma]^2 \[Xi]^3 \[Sigma] +
4 \[Gamma]^3 \[Xi]^3 \[Sigma] +
8 \[Beta] \[Gamma] \[Nu] \[Xi]^3 \[Sigma] +
21 \[Gamma]^2 \[Nu] \[Xi]^3 \[Sigma] +
6 \[Beta] \[Nu]^2 \[Xi]^3 \[Sigma] +
33 \[Gamma] \[Nu]^2 \[Xi]^3 \[Sigma] +
16 \[Nu]^3 \[Xi]^3 \[Sigma] + \[Beta] \[Gamma] \[Xi]^4 \[Sigma] +
2 \[Gamma]^2 \[Xi]^4 \[Sigma] + \[Beta] \[Nu] \[Xi]^4 \[Sigma] +
5 \[Gamma] \[Nu] \[Xi]^4 \[Sigma] +
3 \[Nu]^2 \[Xi]^4 \[Sigma] + \[Beta] \[Gamma]^3 \[Nu] \[Sigma]^2 + \
\[Beta]^2 \[Gamma] \[Nu]^2 \[Sigma]^2 +
5 \[Beta] \[Gamma]^2 \[Nu]^2 \[Sigma]^2 +
3 \[Gamma]^3 \[Nu]^2 \[Sigma]^2 + 2 \[Beta]^2 \[Nu]^3 \[Sigma]^2 +
10 \[Beta] \[Gamma] \[Nu]^3 \[Sigma]^2 +
9 \[Gamma]^2 \[Nu]^3 \[Sigma]^2 + 6 \[Beta] \[Nu]^4 \[Sigma]^2 +
9 \[Gamma] \[Nu]^4 \[Sigma]^2 +
3 \[Nu]^5 \[Sigma]^2 + \[Beta] \[Gamma]^3 \[Xi] \[Sigma]^2 +
2 \[Beta]^2 \[Gamma] \[Nu] \[Xi] \[Sigma]^2 +
7 \[Beta] \[Gamma]^2 \[Nu] \[Xi] \[Sigma]^2 +
7 \[Gamma]^3 \[Nu] \[Xi] \[Sigma]^2 +
5 \[Beta]^2 \[Nu]^2 \[Xi] \[Sigma]^2 +
21 \[Beta] \[Gamma] \[Nu]^2 \[Xi] \[Sigma]^2 +
29 \[Gamma]^2 \[Nu]^2 \[Xi] \[Sigma]^2 +
18 \[Beta] \[Nu]^3 \[Xi] \[Sigma]^2 +
40 \[Gamma] \[Nu]^3 \[Xi] \[Sigma]^2 +
18 \[Nu]^4 \[Xi] \[Sigma]^2 + \[Beta]^2 \[Gamma] \[Xi]^2 \[Sigma]^2 \
+ 2 \[Beta] \[Gamma]^2 \[Xi]^2 \[Sigma]^2 +
4 \[Gamma]^3 \[Xi]^2 \[Sigma]^2 +
4 \[Beta]^2 \[Nu] \[Xi]^2 \[Sigma]^2 +
14 \[Beta] \[Gamma] \[Nu] \[Xi]^2 \[Sigma]^2 +
26 \[Gamma]^2 \[Nu] \[Xi]^2 \[Sigma]^2 +
19 \[Beta] \[Nu]^2 \[Xi]^2 \[Sigma]^2 +
50 \[Gamma] \[Nu]^2 \[Xi]^2 \[Sigma]^2 +
30 \[Nu]^3 \[Xi]^2 \[Sigma]^2 + \[Beta]^2 \[Xi]^3 \[Sigma]^2 +
3 \[Beta] \[Gamma] \[Xi]^3 \[Sigma]^2 +
6 \[Gamma]^2 \[Xi]^3 \[Sigma]^2 +
8 \[Beta] \[Nu] \[Xi]^3 \[Sigma]^2 +
21 \[Gamma] \[Nu] \[Xi]^3 \[Sigma]^2 +
18 \[Nu]^2 \[Xi]^3 \[Sigma]^2 + \[Beta] \[Xi]^4 \[Sigma]^2 +
2 \[Gamma] \[Xi]^4 \[Sigma]^2 +
3 \[Nu] \[Xi]^4 \[Sigma]^2 + \[Beta] \[Gamma]^2 \[Nu] \[Sigma]^3 + \
\[Gamma]^3 \[Nu] \[Sigma]^3 + \[Beta]^2 \[Nu]^2 \[Sigma]^3 +
5 \[Beta] \[Gamma] \[Nu]^2 \[Sigma]^3 +
3 \[Gamma]^2 \[Nu]^2 \[Sigma]^3 + 4 \[Beta] \[Nu]^3 \[Sigma]^3 +
3 \[Gamma] \[Nu]^3 \[Sigma]^3 + \[Nu]^4 \[Sigma]^3 + \[Beta] \
\[Gamma]^2 \[Xi] \[Sigma]^3 + \[Gamma]^3 \[Xi] \[Sigma]^3 +
2 \[Beta]^2 \[Nu] \[Xi] \[Sigma]^3 +
7 \[Beta] \[Gamma] \[Nu] \[Xi] \[Sigma]^3 +
7 \[Gamma]^2 \[Nu] \[Xi] \[Sigma]^3 +
10 \[Beta] \[Nu]^2 \[Xi] \[Sigma]^3 +
14 \[Gamma] \[Nu]^2 \[Xi] \[Sigma]^3 +
8 \[Nu]^3 \[Xi] \[Sigma]^3 + \[Beta]^2 \[Xi]^2 \[Sigma]^3 +
2 \[Beta] \[Gamma] \[Xi]^2 \[Sigma]^3 +
4 \[Gamma]^2 \[Xi]^2 \[Sigma]^3 +
8 \[Beta] \[Nu] \[Xi]^2 \[Sigma]^3 +
15 \[Gamma] \[Nu] \[Xi]^2 \[Sigma]^3 +
14 \[Nu]^2 \[Xi]^2 \[Sigma]^3 + 2 \[Beta] \[Xi]^3 \[Sigma]^3 +
4 \[Gamma] \[Xi]^3 \[Sigma]^3 +
8 \[Nu] \[Xi]^3 \[Sigma]^3 + \[Xi]^4 \[Sigma]^3 + \[Beta] \[Gamma] \
\[Nu] \[Sigma]^4 + \[Beta] \[Nu]^2 \[Sigma]^4 + \[Beta] \[Gamma] \
\[Xi] \[Sigma]^4 +
2 \[Beta] \[Nu] \[Xi] \[Sigma]^4 + \[Gamma] \[Nu] \[Xi] \[Sigma]^4 + \
\[Nu]^2 \[Xi] \[Sigma]^4 + \[Beta] \[Xi]^2 \[Sigma]^4 + \[Gamma] \
\[Xi]^2 \[Sigma]^4 + 2 \[Nu] \[Xi]^2 \[Sigma]^4 + \[Xi]^3 \[Sigma]^4 >0
`

Since all entries are positive on the left hand side, the third condition is satisfied.

Hence by the Routh-Hurwitz condition, all eigenvalues have negative real parts and we conclude that the endemic equilibrium is locally stable.