Symmetrize eigenvectors of degenerate (repeated) eigenvalue

I think I solved it: If matrices $A$ and $S$ commute, i.e. $AS = SA$, then I eigenspaces of one matrix are invariant under the other matrix, because if $Ax = \lambda x$, then $A Sx = SAx = S\lambda x$, so if $x$ is eigenvector of $A$ with eigenvalue $\lambda$, then so is $S\lambda x$.

This in mind, I choose one of the matrices and diagonalize it, say I choose $S$. Then I get $S = T D T^\top$ with $T$ being unitary and $D$ being diagonal. If I then use $T$ to transform $A$, I get $A' = T A T^\top$ is block diagonal, because of the invariance mentioned earlier.

I can then diagonalize each block of $A$ on its own, because that way, the resulting eigenvectors of $A$ will remain eigenvectors of $S$. That way, I have simultaneously diagonalized $A$. The nice thing is that this can help make the resulting eigenvectors unique, because now eigenvectors are characterized not only by their $A$-matrix eigenvalue but also by their $S$-matrix eigenvalue.

If that doesn't suffice, I could search for more symmetries that all simultaneously commute with each other. An example from quantum mechanics is the hydrogen atom, where one find that the Hamiltonian commutes with the modulus of the angular momentum and its z-component, which ultimately gives rise to the three quantum numbers n, l and m.

The solution of Lagerbaer has been tested and i present the results obtained by sympy :

    (A)
    [0  1  1]
    [       ]
    [1  0  1]
    [       ]
    [1  1  0]
    (S)
    [0  1  0]
    [       ]
    [0  0  1]
    [       ]
    [1  0  0]
    eigenvals of A
    [-1  0   0]
    [         ]
    [0   -1  0]
    [         ]
    [0   0   2]
    eigenvects of A
    [-0.707106781186548  -0.707106781186548  0.577350269189626]
    [                                                         ]
    [0.707106781186548           0           0.577350269189626]
    [                                                         ]
    [        0           0.707106781186548   0.577350269189626]
    T.H.Tadjoint
    [-1.0 + 7.46207733321766e-27*I  0.e-137 + 0.e-140*I   -0.e-134 + 0.e-138*I]
    [                                                                         ]
    [     0.e-137 + 0.e-140*I         -1.0 + 0.e-26*I     -0.e-134 + 0.e-138*I]
    [                                                                         ]
    [    -0.e-126 + 0.e-125*I       -0.e-126 + 0.e-125*I          2.0         ]
    eigenvals after TATadj
    [-1.0   0                 0              ]
    [                                        ]
    [ 0    2.0                0              ]
    [                                        ]
    [ 0     0   -1.0 + 7.46207733321766e-27*I]
    eigenvects after TATadj
    [ 0    0   1.0]
    [             ]
    [1.0   0    0 ]
    [             ]
    [ 0   1.0   0 ]

it seems that the method might work