Every way to compute eigenvalues numerically
Power iteration's family, one eigenvalue at a time, vs. the QR algorithm's all-at-once sweep -- and why the characteristic polynomial itself is a trap.
Never form the characteristic polynomial by hand for n>3 -- its coefficients are numerically unstable functions of A's own entries. If you only need the DOMINANT eigenvalue (and its eigenvector), power iteration is the cheapest option, accelerated by the Rayleigh quotient; the inverse power method (shifted near a guess) targets a SPECIFIC eigenvalue instead of the dominant one, and deflation lets you peel off one eigenvalue at a time to reach the next-largest. When you need EVERY eigenvalue, the QR algorithm (after first reducing to Hessenberg form for speed) is the standard, general-purpose method.
Why not the characteristic polynomial
Power iteration and its variants -- one eigenvalue at a time
The single eigenvalue power iteration converges to -- requires it to be strictly larger in magnitude than every other.
Let A be an n\times n diagonalizable matrix with eigenvalues \lambda_1,\lambda_2,\ldots,\lambda_n satisfying \lvert\lambda_1\rvert>\lvert\lambda_2\rvert\ge\lvert\lambda_3\rvert\ge\cdots\ge\lvert\lambda_n\rvert, a single, strictly dominant eigenvalue, with a corresponding eigenbasis \mathbf v_1,\ldots,\mathbf v_n. Suppose a starting vector \mathbf x_0 has a nonzero component along \mathbf v_1. Then the normalized power-iteration sequence \mathbf x_{k+1}=A\mathbf x_k/\lVert A\mathbf x_k\rVert converges, up to sign, to a unit eigenvector for \lambda_1. The Rayleigh quotient \mathbf x_k^{\mathsf T}A\mathbf x_k/\mathbf x_k^{\mathsf T}\mathbf x_k converges to \lambda_1 itself. Both converge at a rate governed by \lvert\lambda_2/\lambda_1\rvert: the error after k steps is O\!\left(\lvert\lambda_2/\lambda_1\rvert^{k}\right).
Show proof
Write \mathbf x_0 in the eigenbasis: \mathbf x_0=c_1\mathbf v_1+c_2\mathbf v_2+\cdots+c_n\mathbf v_n with c_1\ne0 by hypothesis. A\mathbf v_i=\lambda_i\mathbf v_i, so applying A repeatedly gives A^k\mathbf x_0=c_1\lambda_1^k\mathbf v_1+c_2\lambda_2^k\mathbf v_2+\cdots+c_n\lambda_n^k\mathbf v_n =\lambda_1^k\left(c_1\mathbf v_1+\sum_{i=2}^n c_i\left(\frac{\lambda_i}{\lambda_1}\right)^k\mathbf v_i\right). Every ratio \lvert\lambda_i/\lambda_1\rvert<1 for i\ge2. That is exactly what "strictly dominant" means. So (\lambda_i/\lambda_1)^k\to0 as k\to\infty. The term in parentheses converges to c_1\mathbf v_1. The leading scalar factor \lambda_1^k changes only the LENGTH of A^k\mathbf x_0, never its direction. Normalizing at every step removes that factor. (That means dividing by \lVert A^k\mathbf x_0\rVert at each step, exactly what the algorithm does.) This leaves the direction converging to \mathbf v_1 itself, up to an overall sign. (\lambda_1^k can be negative when \lambda_1<0 and k is odd.) The slowest-decaying error term is the one with the largest surviving ratio, \lvert\lambda_2/ \lambda_1\rvert. That ratio governs the overall rate.
The normalized iterate \mathbf x_k converges to a unit eigenvector for \lambda_1. The Rayleigh quotient is continuous in \mathbf x, and exactly equals \lambda_1 when evaluated AT an exact unit eigenvector for \lambda_1 (\mathbf v_1^{\mathsf T}A\mathbf v_1/\mathbf v_1^{\mathsf T}\mathbf v_1=\mathbf v_1^{\mathsf T}(\lambda_1\mathbf v_1)/1=\lambda_1). So continuity carries the vector's own convergence over to the Rayleigh quotient's convergence to \lambda_1, at the same \lvert\lambda_2/\lambda_1\rvert^k rate.
A far more accurate eigenvalue estimate than the raw iterate, from the same iteration.
Power iteration on A^{-1} instead of A -- converges to the SMALLEST eigenvalue.
Shift by a guess near a specific eigenvalue to target that one instead of the smallest.
Term
deflation
Peel off the eigenvalue you've already found so the next power iteration converges to the next-largest.
The QR algorithm -- every eigenvalue at once
The reflection that zeroes out a whole column below the diagonal in one shot -- the workhorse behind Hessenberg reduction.
Almost triangular (one nonzero sub-diagonal) -- reduces every later QR step's own cost from O(n^3) to O(n^2).
Repeatedly factor A_k=Q_kR_k and reassemble as A_{k+1}=R_kQ_k -- converges toward a triangular matrix with the eigenvalues on the diagonal.