Imagine that I have a diagonalizable matrix $A$ (complex-symmetric, btw, but it does not matter for this problem). Some of it's eigenvalues have algebraic, and geometric multiplicities equal to $1$, which is good, but others have geometric multiplicities of $2$ (doubly repeated eigenvalues). This basically means that any linear combination of such eigenvectors is also an eigenvector. For instance, if $\lambda_1 = \lambda_2$ with corresponding $\mathbf{v_1}$, and $\mathbf{v_2}$, then $c_1 \mathbf{v_1} + c_2 \mathbf{v_2}$ is also a valid eigenvector.
Now consider diagonalizable matrix $B$ that is a matrix representation of a symmetry of the system described by the initial matrix $A$, which means that $\left[ A, B \right] = 0$ (they commute), thus, they are simultaneously diagonalizable.
However, due to the degree of freedom mentioned in the beginning, if I will ask a numerical algorithm to return me a set of eigenvectors for $A$ (like function "[V, D] = eig(A)" in MATLAB), they will not necessarily diagonalize $B$ (and vice versa) for that obvious reason. The question is: how I can undergo a procedure of simultaneous diagonalization of $A$, and $B$?
P.S. the reason I want to do this is that the numerical algorithm due to its nature will return me SOME, but not necessarily the nicest and simplest, form of eigenvectors of $A$, but the knowledge of the symmetry $B$ will help me to obtain the simplest, and most convenient form of eigenvectors, hence the question.
All you need is to diagonalise the restriction of $B$ on each eigenspace of $A$. More generally, suppose that $A=PDP^{-1}$ is a diagonalisation. Then one can rearrange the columns of the $P$ so that $$ P^{-1}AP=\pmatrix{\lambda_1 I_{n_1}\\ &\ddots\\ &&\lambda_k I_{n_k}\\ &&&\lambda_{k+1}\\ &&&&\ddots\\ &&&&&\lambda_m} $$ where $\lambda_1,\ldots,\lambda_m$ are distinct and $n_1,\ldots,n_k>1$ (in your case, $n_1=\cdots=n_k=2$). Since $A$ and $B$ commute, $P^{-1}BP$ must be in the form of $$ \pmatrix{B_1\\ &\ddots\\ &&B_k\\ &&&\mu_{k+1}\\ &&&&\ddots\\ &&&&&\mu_m} $$ where each $B_i$ is an $n_i\times n_i$ diagonalisable matrix. Now obtain a diagonalisation $Q_i^{-1}B_iQ_i=D_i$ for each $i$. Then $M^{-1}AM$ and $M^{-1}BM$ are diagonal matrices where $$ M=P\pmatrix{Q_1\\ &\ddots\\ &&Q_k\\ &&&1\\ &&&&\ddots\\ &&&&&1}. $$
Another way to achieve simultaneous diagonalisation is to diagonalise $A+tB$ for an appropriate $t$. If $t$ is chosen at random, then with probability $1$ a diagonalisation of $A+tB$ also diagonalises both $A$ and $B$. So, you use an exclusion method to keep sampling $t$ and diagonalising $A+tB$ until both $A$ and $B$ are simultaneously diagonalised.
If you don't want to pick $t$ at random, you may instead precompute the eigenvalues of $A$ and $B$ and choose any $$ t\not\in\left\{\frac{\lambda-\lambda'}{\mu-\mu'}:\ \lambda,\lambda'\in\operatorname{spectrum}(A), \,\mu,\mu'\in\operatorname{spectrum}(B),\,\mu\ne\mu'\right\}. $$ This condition ensures that no two restrictions of $A+tB$ on two different eigenspaces of $A$ share the same eigenvalue.