Solution (source code)

= Solution

Expanding the <determinant> gives $\det(\lambda I-M)=(\lambda-1)[(\lambda-1)^2-25]$. Solving each <eigenvalue equation> yields
$$
\boxed{\begin{array}{c|c}
\lambda&\text{one eigenvector}\\\hline
6&(4,5,3)^T\\
-4&(4,-5,3)^T\\
1&(3,0,-4)^T
\end{array}.}
$$
Their pairwise <inner products> are $16-25+9=0$, $12-12=0$ and $12-12=0$. Their squared <Euclidean norms> are $50,50,25$, respectively, so division by $\sqrt{50},\sqrt{50},5$ gives an <orthonormal basis>. This also follows from the <spectral theorem for real symmetric matrices>: symmetry makes <eigenvectors> for distinct <eigenvalues> orthogonal.