dS:L-m*S-B*S*(I1+I2);
dI1:B*S*(I1+I2)-m*I1-g*I1;
dR:g*I1-m*R-v*R+r*I2;
dI2:v*R-r*I2-m*I2;

eq:solve([dS,dI1,dR,dI2],[S,I1,R,I2]);

Jac:jacobian([dS,dI1,dR,dI2],[S,I1,R,I2]);

Jac0:Jac,S=L/m,I1=0,I2=0,R=0;

cp:expand(charpoly(Jac0,x));

a1:coeff(cp,x^3);

a2:(coeff(cp,x^2));

a3:coeff(cp,x);

a4:cp,x=0;

/* Routh-Hurwitz stability criteria for n=4   */;
a1>0;
a3>0;
a4>0;
a1*a2*a3 > a3^2+a1^2*a4;