/* LORENZ EQUATIONS (RHS) */
L:[sigma*(y-x),r*x-y-x*z,x*y-b*z];

solve(%,[x,y,z]);

ratsimp(%);

ratsimp(substitute([x=sqrt(b*r-b),y=sqrt(b*r-b),z=r-1],L));

J:matrix([diff(L[1],x), diff(L[1],y), diff(L[1],z)],[diff(L[2],x), diff(L[2],y), diff(L[2],z)],[diff(L[3],x), diff(L[3],y), diff(L[3],z)]);

J:substitute([x=sqrt(b*r-b),y=sqrt(b*r-b),z=r-1],J);

ratsimp(eigenvalues(J));

load(facexp);

px:collectterms(expand(-charpoly(J,x)),x);

/* coefficients for R-H criteria */
c3:substitute(x=0,px);

c2:coeff(px,x);

c1:coeff(px,x^2);