/* Linear 2D example with a nonlinear perturbation in the 1st equation */
f: a*x+b*y;
g: c*x+d*y;

/* Solving for f=0, g=0 we get the Eq. Pt. (0,0). */
eq0: solve([f,g],[x,y]);

/*  The Jacobian at 0,0 is... */
J: jacobian([f,g],[x,y]);
J0:J,eq0;

/* Characteristic Polynomial for the Routh-Hurwitz criteria */
p: determinant(J0-matrix([L,0],[0,L]));

p0:expand(p);

c1:coeff(p0,L);

c2:p0, L=0;

/* Stability Criteria: c1>0, c2>0 */
c1>0;
c2>0;

/*NOTE: This result applies for any model generated by adding a nonlinear term that preserves (0,0) as an equilibrium point. */
f2: a*x+b*y+c*y^3;

J2: jacobian([f2,g],[x,y]);

J2, x=0, y=0; /* evaluate at (0,0) */;

/* same jacobian, same eigenvalues, etc., so see above! */ 