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

dN;

solve(dN,P);

eqs:solve([dN,dP],[N,P]);

/* When is that last one positive? */
eq2:eqs[4];

eq2[1]>0;

/* positive if ac > mu */;

eq2[2]>0;

factor(num(rhs(eq2[2])))/(c*k*r)>0;

/* acK > mu*(k+K) */;

/* Now let's find our Jacobian... */
Jac:(jacobian([dN,dP],[N,P]));

ratsimp(Jac);

Jac;

Jac, [N=0,P=0];

Jac, eqs[2];

Jac, eq2;

J2:ratsimp(%);

eigenvalues(J2);

Tr:J2[1,1];

Det:determinant(J2);

/* Stability criteria: Tr < 0 , Det > 0 */
Tr<0;
ratsimp(-Det*a*c*K/mu/r)-K*(mu-a*c)<-K*(mu-a*c);

/* Since we're assuming N>0 aka ac-mu > 0 */

Tr<0;
(ratsimp(-Det*a*c*K/mu/r)-K*(mu-a*c))/(K*(a*c-mu))<1;

ratsimp((ratsimp(-(factor(Tr)*a*c*K/mu/r)*(mu-a*c))+k*mu)/((a*c-mu)*K))<k*mu/((a*c-mu)*K);

fsff

/* So we have ... */
1-(a*c*k)/((a*c-mu)*K)<(k*mu)/((a*c-mu)*K);

/* So our 1st staility condition, assuming N>0 is */
N<K;
/* Second... ... */
1-(a*c)*N/(mu*K)<N/K;

1/(1+(a*c)/(mu))<N/K;

Tr;