Numerical method for the Eigenvectors of a $3\times3$ symmetric matrix I am looking for an ad-hoc algorithm to quickly find the Eigenvectors of a symmetric $3\times3$ matrix. I already have a solution for the Eigenvalues by direct resolution of the characteristic equation.
I want a solution that is efficient (presumably non-iterative), yet is able to handle corner cases such as equal Eigenvalues. An option could be to unroll a general algorithm to the case $3\times3$, but this can get complicated.
As the Eigenvalues are available Gaussian elimination is possible, but I don't see how it could handle equal Eigenvalues.
Any suggestion ?

Update:
Vectorially, when $\lambda$ is an Eigenvalue the three column vectors of $M-\lambda I$ are coplanar. Then the coefficients of the linear combination of those three vectors correspond to the vector normal to the plane.
Numerically, one can form the three cross-products of the vectors in pairs and keep the product with the largest norm. This involves the computation of nine $2\times2$ minors and three dot squares ($15$ additions and $27$ multiplications), and must be repeated for the three Eigenvalues.
When two Eigenvalues are equal, the vectors are parallel so that the cross-products cancel.

Second update:
A dream article: "A Robust Eigensolver for 3 × 3 Symmetric Matrices
David Eberly". https://www.geometrictools.com/Documentation/RobustEigenSymmetric3x3.pdf
 A: The problem you are talking about is numerically unstable: at the identity matrix, any orthogonal basis could be your answer, but by arbitrarily small perturbations of the identity matrix you could force any one particular orthogonal basis to be the unique (up to scalar factors $-1$) answer. Therefore your goal of having a single (continuous) formula, parametrised by the eigenvalue, describing the eigenvectors is unattainable.
Also it is not clear to me how you would go about having a formula in the case of multiple eigenvalues; do you imagine having a separate formula for the first (chosen) eigenvalue at$~\lambda$, then another formula for the second one, valid only if $\lambda$ is a multiple eigenvalue? In any case the stability argument shown no continuous formula for the first eigenvalue can exist.
A: Gaussian elimination will give you a non-orthogonal basis for an eigenspace even in the case where there are equal eigenvalues.  You can then orthogonalize this basis by Gram-Schmidt.  You will have to pick some tolerances for deciding when two or three eigenvalues are equal.  
