$Post:=If[MatrixQ[#],MatrixForm[#],#]&
<<LinearAlgebra`Master`

Exercise 3.3.1

M=Array[Max,{3,3}]

Id[3]=IdentityMatrix[3];
p=Det[x Id[3]-M]

lambda=Solve[p==0,x][[1,1,2]]

q=p/.{x->r};
CM=r Id[3]-M

CM[[1]]=CM[[1]]/CM[[1,1]]; CM

CM[[2]]=CM[[2]]+2CM[[1]];
CM[[3]]=CM[[3]]+3CM[[1]];
CM=CM//Together

CM[[2]]=CM[[2]]/CM[[2,2]];
CM=CM//Together

CM[[3]]=CM[[3]]-CM[[3,2]]CM[[2]];
CM[[1]]=CM[[1]]-CM[[1,2]]CM[[2]];
CM=CM//Together; CM//TableForm

CM[[3,3]]=0; CM

NullSpace[CM]

lambda=Solve[p==0,x];

Do[r[i]=lambda[[i,1,2]],{i,1,3}]

Do[ES[i]=%%%/.{r->r[i]},{i,1,3}]

P=Transpose[Join[ES[1],ES[2],ES[3]]];

NullSpace[r Id[3]-M,ZeroTest->(PolynomialMod[#,q]==0&)] \
         //Together

In order not to repeat parts of the exercise, the equivalence
of the above result with NullSpace[CM] will be shown at the
end of the next exercise.

NS=%

Exercise 3.3.2

P=N[P];
Unprotect[D];
D=Inverse[P].M.P;
D1[n_]:=Table[If[Abs[Re[D[[i,j]]]]<=0.000001,                 \
                 I Im[D[[i,j]]],D[[i,j]]],{i,n},{j,n}]
Diag[n_]:=Table[If[Abs[Im[D1[n][[i,j]]]]<=0.000001,           \
                   Re[D1[n][[i,j]]],D1[n][[i,j]]],{i,n},{j,n}]
Diag[3]

Equivalence from previous exercise:

Do[ES[i]=NS/.{r->r[i]},{i,1,3}]
P=Transpose[Join[ES[1],ES[2],ES[3]]];
P=N[P];
D=Inverse[P].M.P;
Diag[3]

Exercise 3.3.3

M=Array[Max,{4,4}]

x IdentityMatrix[4]-M

CM=Join[RowReduce[Array[%[[#]]&,3]],{{0,0,0,0}}];
CM=CM//Together; CM//TableForm

Clear[r];
p=Det[r IdentityMatrix[4]-M]

Solve[p==0,r];
r=r/.%;
P=Transpose[Table[NullSpace[CM][[1]]/.{x->r[[i]]//N},{i,4}]//N]//Re

Although it is not specified in the exercises, we will use
floating point representations in this chapter, since exact
calculations here would overflow the capacity of personal
computers at the present time. 

D=Inverse[P].M.P;
Diag[4]

Exercise 3.3.4

M=Array[Min,{4,4}]

x IdentityMatrix[4]-M

CM=Join[RowReduce[Array[%[[#]]&,3]],{{0,0,0,0}}];
CM=CM//Together; CM//TableForm

Clear[r];
p=Det[r IdentityMatrix[4]-M]

Solve[p==0,r];
r=r/.%;
P=Transpose[Table[NullSpace[CM][[1]]/.{x->r[[i]]//N},{i,4}]//N]//Re

Inverse[P].M.P;
Diag[4]

Exercise 3.3.5

M=Array[Max,{5,5}]

x IdentityMatrix[5]-M

CM=Join[RowReduce[Array[%[[#]]&,4]],{{0,0,0,0,0}}];
CM=CM//Together; CM//TableForm

Clear[r]
p=Det[r IdentityMatrix[5]-M]

NSolve[p==0,r];
r=r/.%;
P=Transpose[Table[NullSpace[CM][[1]]/.{x->r[[i]]//N},{i,5}]//N]//Re;
P//Together//TableForm

D=Inverse[P].M.P;
Diag[5]

Exercise 3.3.6

M=Table[If[i>3,If[j>3,Min[i,j],0],If[j>3,0,Max[i,j]]], \
                                           {i,6},{j,6}]

Trying to solve for the 6*6 matrix as a whole, we fail. 
However, we notice that by finding diagonalizing matrices 
for the upper left and the lower right corner of M we can 
then construct the diagonalizing matrix of M by the same 
pattern M was constructed.

M1=SubMatrix[M,{1,1},{3,3}]

M2=SubMatrix[M,{4,4},{3,3}]

x Id[3]-M1

CM=Join[RowReduce[Array[%[[#]]&,2]],{{0,0,0}}];
CM=CM//Together; CM//TableForm

Clear[r]
p=Det[r Id[3]-M1]

Solve[p==0,r];
r=r/.%

P=Transpose[Table[NullSpace[CM][[1]]/.{x->r[[i]]//N},{i,3}]//N]//Re

D=Inverse[P].M1.P;
Diag[3]

P1=P

x Id[3]-M2

CM=Join[RowReduce[Array[%[[#]]&,2]],{{0,0,0}}];
CM=CM//Together; CM//TableForm

Clear[r]
p=Det[r Id[3]-M2]

Solve[p==0,r];
r=r/.%

P=Transpose[Table[NullSpace[CM][[1]]/.{x->r[[i]]//N},{i,3}]//N]//Re

D=Inverse[P].M2.P;
Diag[3]

P2=P

ZeroMatrix=Array[0&,{3,3}]

P=BlockMatrix[{{P1,ZeroMatrix},{ZeroMatrix,P2}}]

D=Inverse[P].M.P;
Diag[6]//TableForm