I have a Hermitian Matrix $H$ that evolves very slightly over ~50 iterations until convergence of a Schrödinger/Poisson system is achieved. The eigenvalues I'm interested in are $\lambda_1$ and the degenerate $\lambda_2 = \lambda_3$. The eigenvectors $v_n$ are all orthogonal. I'm interested in the evolution of the Eigenvectors for which I have very accurate starting guesses $v_{n}$. I'm calculating $\lambda_n$ using Rayleigh-coefficients and plug those into the inverse iteration $(H-\lambda_n I )^{-1} v_{n} = v_{n}'$.
My problem is, that there is no way of knowing which of the degenerate Eigenvectors I get when plugging in $v_2$ and $v_3$. I've tried using Lanczos-algorithm to calculate the first three eigenvectors and eigenvalues and then sort them manually using overlap integrals, however this approach runs into a myriad of problems in the long run which is why I've resorted to calculating the Eigenvectors myself.
My question now is: Is there any way to update degenerate eigenvectors e.g. $v_2$ to $v_2'$ and $v_3$ to $v_3'$ without manually sorting them afterwards? I'm more interested in updated eigenvectors similar to the old eigenvector than accuracy of eigenvectors or eigenvalues.
A similar question was asked here, however being an engineer and not a mathematician I was not able to do anything with the suggested approaches How to get the same eigenvectors of a degenerate eigenvalue as the matrix evolves slightly?