A=sym([3 1 1 0;-1 1 0 -1;-1 -1 0 1;-2 -1 -1 1]) % Use symbolic math toolbox
evals = eig(A) % the eigenvalues are 1 1 1 2
% Find an eigenvector for eval 2
rrefAminus2I = rref(A - 2*eye(4) )  % one free variable
vlam2 = [-1; 0; 1; 1]  % a column vector  
% Note that I don't need sym( ) in previous line, but it wouldn't hurt.
shouldBe0 = rrefAminus2I*vlam2
% Now handle the eigenvalue 1, with multiplicity 3.
AminusI = A - eye(4)
rrefAminusI = rref(AminusI) % result tells me delta_1 = 2
rrefAminusI2 = rref(AminusI^2) % result tells me delta_2 = 3 (multiplicity)
v12 = [0; 0; 0; 1]
shouldBe0 = AminusI^2 *v12
v11 = AminusI*v12  % nonzero so I'm OK
% Now choose a second eigenvector LI from v11
v21 = [-1; 2; 0; 1] % obtained from free variables 0, 1 in rrefAmI
shouldBe0 = AminusI*v21
syms t c1 c2 c3 c4;
x(t) = vlam2*exp(2*t)*c1 + exp(t)*(v11*(c2 + c3*t) + v12*c3 + v21*c4)
ODEcheckShouldBe0 = simplify(diff(x, t)- A*x)
% Now comute the JCF of A, which should be diag(J1(2), J2(1), J1(1))
P = [vlam2 v11 v12 v21]
B = P^-1*A*P  % Indeed, B is in JCF