Attempt to create a perpendicular matrix P for M.

$>$
load(eigen)$$>$ M: matrix([+1,0, 2], [0, -2, +1], [2, +1, -1]); $\left(\begin{array}{ccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 2\hfill \\ \hfill 0\hfill & \hfill -2\hfill & \hfill 1\hfill \\ \hfill 2\hfill & \hfill 1\hfill & \hfill -1\hfill \end{array}\right)\text{}$ $>$ PL: uniteigenvectors(M); $\left[\left[\left[-\frac{\sqrt{13}-1}{2},\frac{\sqrt{13}+1}{2},-3\right],\left[1,1,1\right]\right],\left[\left[\left[\frac{2}{\sqrt{2\sqrt{13}+13}},-\frac{\sqrt{13}+3}{2\sqrt{2\sqrt{13}+13}},-\frac{\sqrt{13}+1}{2\sqrt{2\sqrt{13}+13}}\right]\right],\left[\left[\frac{2}{\sqrt{13-2\sqrt{13}}},\frac{\sqrt{13}-3}{2\sqrt{13-2\sqrt{13}}},\frac{\sqrt{13}-1}{2\sqrt{13-2\sqrt{13}}}\right]\right],\left[\left[\frac{1}{3},\frac{2}{3},-\frac{2}{3}\right]\right]\right]\right]\text{}$ $>$ PLm: apply(’matrix, PL); $\left(\begin{array}{c}\hfill \left[\frac{2}{\sqrt{2\sqrt{13}+13}},-\frac{\sqrt{13}+3}{2\sqrt{2\sqrt{13}+13}},-\frac{\sqrt{13}+1}{2\sqrt{2\sqrt{13}+13}}\right]\hfill \\ \hfill \left[\frac{2}{\sqrt{13-2\sqrt{13}}},\frac{\sqrt{13}-3}{2\sqrt{13-2\sqrt{13}}},\frac{\sqrt{13}-1}{2\sqrt{13-2\sqrt{13}}}\right]\hfill \\ \hfill \left[\frac{1}{3},\frac{2}{3},-\frac{2}{3}\right]\hfill \end{array}\right)\text{}$ The problem with PL is, that it doesn’t have a list of eigenvectors as its second output element. It has a list, of lists, of eigenvectors, sorted at top-level per eigenvalue. There just happens to be one eigenvector per eigenvalue in this case. This is why, to convert that to a matrix, results in a column vector of lists. In turn, trying to transpose that, results in a row vector, of lists. There’s nothing wrong in the way the matrix is being built. It’s just useless to me, how the uniteigenvectors() function outputs its eigenvectors. $>$ P: ﬂoat(transpose(PLm)); $\left(\begin{array}{ccc}\hfill \left[0.4448719185337484,-0.7346560672221786,-0.5122201079553044\right]\hfill & \hfill \left[0.8312507834516553,0.1258412430373975,0.5414666347632251\right]\hfill & \hfill \left[0.3333333333333333,0.6666666666666666,-0.6666666666666666\right]\hfill \end{array}\right)\text{}$ $>$ ReduceDepth(L) := block ( Output : [], for i in L do ( Output : append(i, Output) ), return(Output) )$

$>$
P: ﬂoat(transpose(apply(’matrix, ReduceDepth(PL))));
$\left(\begin{array}{ccc}\hfill 0.3333333333333333\hfill & \hfill 0.8312507834516553\hfill & \hfill 0.4448719185337484\hfill \\ \hfill 0.6666666666666666\hfill & \hfill 0.1258412430373975\hfill & \hfill -0.7346560672221786\hfill \\ \hfill -0.6666666666666666\hfill & \hfill 0.5414666347632251\hfill & \hfill -0.5122201079553044\hfill \end{array}\right)\text{}$

The order of columns has been reversed, but
the result will be equivalent.

$>$
D: transpose(P).M.P;
$\left(\begin{array}{ccc}\hfill -3.0\hfill & \hfill 0.0\hfill & \hfill -1.110223024625156\cdot {10}^{-16}\hfill \\ \hfill 0.0\hfill & \hfill 2.302775637731994\hfill & \hfill -2.220446049250313\cdot {10}^{-16}\hfill \\ \hfill 0.0\hfill & \hfill 0.0\hfill & \hfill -1.302775637731994\hfill \end{array}\right)\text{}$