J: matrix(
 [a,b], 
 [c,d]
);

p:determinant(J-matrix([L,0],[0,L]));

expand(p);

eigval:solve(%,L);

p:L^2-TrJ*L+DetJ;

solve(p,L);

v:matrix([L-d],[c]);

pi:[a=1/10,b=-1,c=1/2,d=-4];
Li:ratsimp(substitute(pi,eigval));
ratsimp(substitute([Li[1],pi[3],pi[4]],v));
ratsimp(substitute([Li[2],pi[3],pi[4]],v));

pii:[a=1,b=1,c=-1,d=-2];
Lii:ratsimp(substitute(pii,eigval));
ratsimp(substitute([Lii[1],pii[3],pii[4]],v));
ratsimp(substitute([Lii[2],pii[3],pii[4]],v));

piii:[a=-1,b=1,c=-1,d=-1];
Liii:ratsimp(substitute(piii,eigval));
ratsimp(substitute([Liii[1],piii[3],piii[4]],v));
ratsimp(substitute([Liii[2],piii[3],piii[4]],v));

/* Part c */
Li,numer;
Lii,numer;
Liii,numer;

/* BONUS */
M: matrix(
 [a,b,c], 
 [d,e,f], 
 [g,h,i]
);

cp:determinant(M-matrix([x,0,0],[0,x,0],[0,0,x]));

load(facexp);

/* Rearange */
collectterms(expand(cp),x);

/* a+e+i = Tr(M) */
/* constant term, by definition, is Det(J) */
/* The x coefficient is-1/2(Tr(A)^2-Tr(A^2)) */
expand(-1/2*(mat_trace(M)^2 - mat_trace(M.M)));

ratsimp(expand(-1/2*(mat_trace(M)^2 - mat_trace(M.M))) - coeff(expand(cp),x));

/* Thus we could write ... */
cp:-x^3-(-TrM)*x^2- (-1/2*(TrM^2 - TrMM))*x - DetM;

/* Since cp(x)=0 has the same solutions as -cp(x)=0, rewrite as p=-cp ... */
p: x^3 + (-TrM)*x^2 + (-1/2*(TrM^2 - TrMM))*x + DetM;

/* Now we can lookup a general formula for the roots of a (monic) cubic equation */
/* Writing things in term of the discriminant D ... */
D:  18*(-TrM)*(TrMM-TrM^2)/2*(DetM);
/* ... and we additionally write ... */
D0: TrM^2-3*(-1/2*(TrM^2 - TrMM));
D1: -2*TrM-9*TrM*1/2*(TrM^2 - TrMM)+27*DetM;
C:  (D1+sqrt(-27*D)/2)^(1/3);

/* The roots of our cubic equation (i.e., eigenvalues of M) are ... */
L1=ratsimp(-1/3*(TrM+C+D0/C));
L2=ratsimp(-1/3*(TrM+C*((-1+sqrt(-3))/2)+D0/C/((-1+sqrt(-3))/2)));
L2=ratsimp(-1/3*(TrM+C*((-1-sqrt(-3))/2)+D0/C/((-1-sqrt(-3))/2)));