# Is the QR Algorithm guaranteed to compute eigenvectors?

I'm writing some C++ matrix library for hobby. For computing eigenvalues and eigenvectors, I referred the following "Francis double step QR algorithm":

In particular, page 82 of https://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf

My implementation correctly computes eigenvalues, but fails to compute eigenvectors, except the first eigenvector corresponding to the largest eigenvalue.

For $$\begin{pmatrix} -2 & 2 & 3 & 4 \\ -9 & 7 & 5 & 5 \\ -5 & 2 & 6 & 6 \\ -7 & 2 & 8 & 9 \end{pmatrix}$$, Wolframalpha gives the following result:

$$\lambda_{1} = 12.6009, \lambda_{2} = 4.2635, \lambda_{3} = 2.5983, \lambda_{4} = 0.5373$$ with corresponding eigenvectors

$$v_1 = (0.520545, 0.706921, 0.728858, 1)$$

$$v_2 = (-2.02909, -7.68532, -0.446183, 1)$$

$$v_3 = (8.65707, 12.3871, 3.67796, 1)$$

$$v_4 = (-0.748529, -0.608655, -1.56064, 1)$$

My code outputs the following:

(12.6009,0) {(0.343092,0), (0.465934,0), (0.480393,0), (0.659103,0)}
Av(normalized) : {(0.343092,0), (0.465934,0), (0.480393,0), (0.659103,0)}
(4.26349,0) {(0.101326,0), (0.836911,0), (-0.197859,0), (-0.500164,0)}
Av(normalized) : {(-0.180481,0), (0.234039,0), (-0.48548,0), (-0.822777,0)}
(2.59833,0) {(0.924778,0), (-0.243128,0), (-0.270181,0), (-0.112592,0)}
Av(normalized) : {(-0.203286,0), (-0.674777,0), (-0.41863,0), (-0.572799,0)}
(0.537284,0) {(0.12964,0), (-0.152886,0), (0.810603,0), (-0.55022,0)}
Av(normalized) : {(-0.276698,0), (-0.77434,0), (0.503774,0), (0.264664,0)}


,0 part is the imaginary part.

I collected the column vectors of the $$Q$$ matrix upon convergence of the algorithm.

The eigenvalues are correct, and the first eigenvector $$\mathbf{v}$$ is also correct - being equal to the normalized $$\mathbf{A}\mathbf{v}$$. But other eigenvectors seem just broken. Really weird that only one eigenvector is correct.

Is the paper wrong? Should I use another algorithm for computing eigenvectors? Or just is my implementation broken?