/* TREE COMPETITION MODEL */
dH:rH*H-aH*H^2-bH*S*H;
dS:rS*S-aS*S^2-bS*H*S;

/* Find equilibrium points */
EQ:solve([dH=0,dS=0],[H,S]);

/* Stability? First, find the general Jacobian 
   Entries are just the partials J[i,j]=diff(dxi,xj) 
   Instead of typing all that in, just use jacobian() */
J:jacobian([dH,dS],[H,S]);

/* Evaluate the Jacobian at each EQ point */
substitute(EQ[1],J);

/* eigenvalues are rH>0 and rS>0, so (0,0) always unstable */
eigenvalues(%o5);

substitute(EQ[2],J);

eigenvalues(%o7);
/* As above, ignore the "[1,1]" */;

/* Thus, criteria for the stability of EQ[2] are */
num(%o8[1][1])>0;

/* The other 1-species equilibrium should be similar, by symmetery */
substitute(EQ[3],J);

eigenvalues(%o10);

eigenvalues(substitute(EQ[3],J));

/* Now, the coexistence equilibrium... */
EQ[4];

substitute(EQ[4],J);

/* Note that just jumping into finding eigenvalues gets messy! */
eigenvalues(substitute(EQ[4],J));

/* HERE ARE THREE APPROACHES TO GETTING SIMPLER STABILTY CRITERIA */
/* #1 -- simplify a little */
ratsimp(substitute(EQ[4],J));

/* This rarely helps much... */
eigenvalues(%);

/* #2 -- use "f(x_equilbrium)=0" and leave (Seq,Heq) in J to simplify!  */
/* Assuming H!=0 and S!=0... */
ratsimp(dH/H)=0;
ratsimp(dS/S)=0;

J;

matrix([-bH*S+rH-aH*H -aH*H,-bH*H],[-bS*S,rS-bS*H-aS*S-aS*S]);

/* Thus, we have the slightly nicer looking Jacobian */
J4:matrix([0 -aH*H,-bH*H],[-bS*S,0-aS*S]);

/* Now, without substituting EQ4=(H,S) into J4... */
eigenvalues(J4);

/* Better! But still messy. See formula from problem 1 for comparison. */
ratsimp(%);

/* Alternatively, give up on finding eigenvalues, we can use the Routh-Hurwitz criteria instead! */
/* First, get the characteristic polynomial of J(EQ[4])... */
substitute(EQ[4],J);

/* We can do it in terms of paramter values.. */
ratsimp(%);

/* ... or for now leaving the equilibrium expressions for H and S out of J */
J4:ratsimp(J4);

/* Characteristic Polynomial */
cp:charpoly(J4,x);

/* To get a monic polynomial with properly defined coeffients, get the collectterms() function: */
load(facexp);

cp:collectterms(expand(cp),x);

/* Since it's monic (i.e., the x^2 coefficient is +1) */
/* this is of the form p(x)=x^2+c1*x+c0 */
c0:factor(substitute(x=0,cp));
c1:factor(coeff(cp,x));

/* Therefore, we have stability criteria... */
c0>0;

/* which, if H*S>0 is the same as */
aH*aS - bH*bS > 0;

/* and also the second criterion */
c1>0;

/* which always holds for EQ4 > 0. Recall that */
EQ[4];

substitute(EQ[4][2],c0>0);

/* Therefore, if we assume we only care about the case
   where EQ[4]>0 (H>0, S>0), then c0>0 and c1>0 (the Routh-Hurwitz 
   stability criteria) are equivalent to... */
EQ[4][1]>0;
EQ[4][2]>0;
aH*rS-bS*rH>0;

/* Which together are equivalent to to the conditions
   given in the text on page 119: */
aH*aS-bH*bS>0;
rH*aS-rS*bH>0;
aH*rS-bS*rH>0;