O15-C914
> restart:with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> a:=matrix(3,3,[1,3,-7,2,6,-14,1,3,-7]);

                              [1    3     -7]
                              [             ]
                         a := [2    6    -14]
                              [             ]
                              [1    3     -7]

> eigenvects(a);

                   [0, 3, {[-3, 1, 0], [7, 0, 1]}]

> colspace(a);#e1=[7, 0, 1],e2=[1, 2, 1],e3=[1,0,0]

                             {[1, 2, 1]}

> p:=transpose(matrix(3,3,[[7, 0, 1],[1, 2, 1],[1,0,0]]));

                               [7    1    1]
                               [           ]
                          p := [0    2    0]
                               [           ]
                               [1    1    0]

> t:=evalm(p^(-1)&*a&*p);

                               [0    0    0]
                               [           ]
                          t := [0    0    1]
                               [           ]
                               [0    0    0]

> z:=matrix(3,3):
> sys:=convert(evalm(z&*t-t&*z),set);

 sys := {z[1, 2], z[3, 2], z[2, 2] - z[3, 3], -z[3, 1], -z[3, 2], 0}

> s:=solve(sys,{z[1,2],z[3,2],z[3,1],z[2,2]});

   s := {z[1, 2] = 0, z[3, 1] = 0, z[3, 2] = 0, z[2, 2] = z[3, 3]}

> assign(s):evalm(z);

                   [z[1, 1]       0       z[1, 3]]
                   [                             ]
                   [z[2, 1]    z[3, 3]    z[2, 3]]
                   [                             ]
                   [   0          0       z[3, 3]]

> sys1:=convert(evalm(z&*z-t),set);

  sys1 := {z[1, 1] z[1, 3] + z[1, 3] z[3, 3],

                                                     2         2
        z[2, 1] z[1, 1] + z[3, 3] z[2, 1], 0, z[3, 3] , z[1, 1] ,

        z[2, 1] z[1, 3] + 2 z[3, 3] z[2, 3] - 1}

> s1:=solve(sys1);

                                                      1
  s1 := {z[2, 1] = z[2, 1], z[1, 1] = 0, z[1, 3] = -------,
                                                   z[2, 1]

        z[2, 3] = z[2, 3], z[3, 3] = 0}

> assign(s1):evalm(z);

                   [                         1   ]
                   [   0          0       -------]
                   [                      z[2, 1]]
                   [                             ]
                   [z[2, 1]    z[3, 3]    z[2, 3]]
                   [                             ]
                   [   0          0          0   ]

> evalm(z^2);

                            [0    0    0]
                            [           ]
                            [0    0    1]
                            [           ]
                            [0    0    0]

>