/* Model equations */
dN : r*N*(1-N/K) - a*N*P/(1+a*h*N);
dP : c*a*N*P/(1+a*h*N) - mu*P;

/* substituting n=N/K, p=c*P/K yields */
dn : dN/K,N=n*K,P=p*K*c,a=1/k/K/h;
dp : dP/K/c,N=n*K,P=p*K*c,a=1/k/K/h;

dn:dn,c=b*h;
dp:dp,c=b*h;

dn:expandwrt(dn,K); /* r would also work */;

dp:expandwrt(dp,h,k);

/* Now our scaled equations are */
dn;
dp;

/* We could further rescale time tau=mu*t  */
/* This gives the following new parameters */
/* These differ from the above reparameterization! */
/* Here                                            */
/*    p=P/K/c, n=N/K, g=r/mu, k=1/a/K/h, b=c/h/mu  */
/* */
dn : dN/K/mu,N=n*K,P=p*K*c,a=1/k/K/h,r=g*mu;
dp : dP/K/c/mu,N=n*K,P=p*K*c,a=1/k/K/h;

dn:dn,c=b*h*mu;
dp:dp,c=b*h*mu;

dn:expandwrt(dn,K);
dp:expandwrt(dp,h,k);

/* ********************* */
/* Next, equlibria....   */
/* ********************* */
EQs: solve([dn,dp],[n,p]);

EQ0:EQs[1];

EQn:EQs[2];

/* Ignore that n=-k case if focusing on n,p>0 */
EQnp:EQs[4];

/* This gives us criteria for n>=0,p>=0. */
/* b>1 and... */
rhs(EQs[4][2])>0;

/* if b>1, g>0, then the denom b^2-2b+1 = (b-1)^2 > 0 thus */
factor(num(rhs(EQs[4][2]))/g)>0;

factor(num(rhs(EQs[4][2]))/g)/k>0;

factor(num(rhs(EQs[4][2]))/g)/k + k > k;

/* Note: This is also the requirement for n=k/(b-1) to be in (0,1)! */
/* **************************************************************** */
/* Now we determine stability criteria for each EQ */
J : jacobian([dn,dp],[n,p]);

/* first, the origin */
J0:J,n=0,p=0;

/* Therefore we never have stability unless g=r/mu < 0 */
eigenvalues(%);

/* Next, the predator-free equilibrium */
Jn:J,n=1,p=0;

eigenvalues(Jn); /* these are the diagonal entries */;

/* And the coexistence equilbrium. Here we apply two approaches: */
/* First, we just use the equlibrium values and simiplify... */
Jnp : J,EQnp[1],EQnp[2];

Jnp:ratsimp(Jnp);

/* Next, we further simplify using our equlibrium values */
EQnp;

/* Note that ... */
p=factor(rhs(EQnp[2]));

%,solve(EQnp[1],k);

pofn:factor(ratsimp(%));

/* This should be obvious from dn,dp, and/or from solving dn+dp=0 for p! */
/* Simalarly... */
ratsimp(dp/p+1)=1;

J,EQnp[1],ratsimp;

J,pofn,ratsimp;

J,pofn,ratsimp;

%,EQnp[1],ratsimp;

%,solve(EQnp[1],k),ratsimp;

/* Manually, partial substitution gives... */
J2:matrix([(-g*b*n+g*(1-n))/b,-1],[factor(-((b-1)*g*n+(1-b)*g)/b),0]);

cp:charpoly(J2,L); /* See the Routh-Hurwitz Criteria */;

/* Or alternatively, since dim(J)=2 */
determinant(J2) >0 ;
J2[1,1]+J2[2,2] <0;

/* Since b>1 by n>0, and n<1 by p>0, we have */
/* that the first condition holds. */
/* The second condition yields... */
(J2[1,1]+J2[2,2])/g*b <0, ratsimp;

/* Which rearranges to */
(lhs(%) + n*(b+1))/(b+1) < (rhs(%) + n*(b+1))/(b+1), ratsimp;

/* ********************************************* */
/* Alternatively, use the Routh-Hurwitz Criteria */
cp:expand(cp);

coeff(cp,L);

cp,L=0,factor;