How To Actually Get A Matrix Into Diagonal Form Without Losing Your Mind
I spent way too many hours in grad school wrestling with Jordan normal form because someone forgot to tell me that diagonalization is only the clean case. Here is what actually works when you are handed a 4x4 matrix and told to diagonalize it, plus the gotchas I keep running into even now. A matrix A is diagonalizable if you can find an invertible matrix P and a diagonal matrix D such that A = PDP^(-1). That means the columns of P are eigenvectors of A and the diagonal entries of D are the corresponding eigenvalues. That is the entire definition. Everything else is mechanics. Set the characteristic polynomial to zero. That is det(A - lambda * I) = 0. For a 2x2 or 3x3, this is usually straightforward. For a 4x4 or larger, you are going to need a computer or some very clever structure in the matrix. Do not try to expand a 5x5 determinant by hand unless you enjoy pain.
Once you have the polynomial, factor it completely. Repeated roots are not the end of the world, but they change the game later. If your characteristic polynomial is (lambda - 2)^3 * (lambda + 1), you have eigenvalue 2 with algebraic multiplicity 3 and eigenvalue -1 with algebraic multiplicity 1. The number 2 will need three linearly independent eigenvectors for the matrix to be diagonalizable. If it only gives you two, you are stuck with a Jordan block and the whole diagonalization approach collapses.
Step Two: Find The Eigenvectors For Each Eigenvalue
For each eigenvalue lambda, solve the homogeneous system (A - lambda * I) * v = 0. The null space of that matrix gives you the eigenvectors. Row reduce A - lambda * I. The free variables in the reduced form directly give you the basis vectors for the eigenspace. Here is where most people mess up. They find one eigenvector and assume they are done. You need a full basis for the eigenspace, and the number of vectors in that basis must equal the algebraic multiplicity of the eigenvalue. If the geometric multiplicity is less than the algebraic multiplicity, the matrix is not diagonalizable and you need to move on to other methods. I once worked on a control systems problem where the system matrix looked perfectly reasonable, and I spent twenty minutes convinced I had made an arithmetic error because the eigenvalue 5 had algebraic multiplicity 2 but only yielded one eigenvector. Turned out the matrix genuinely was defective. The workaround was to compute the generalized eigenvector by solving (A - 5I)v2 = v1, where v1 is the ordinary eigenvector. That gave me the Jordan chain I needed, and the matrix exponential followed from the Jordan form instead. It added about ten extra minutes of work but saved me from reporting a wrong answer.
Get the Full Details

Step Three: Assemble P And D
Put all the eigenvectors as columns of P. Put the corresponding eigenvalues on the diagonal of D. The order matters, and you have to keep them matching. If the first column of P is an eigenvector for eigenvalue 3, then D[1,1] must be 3. Mix that up and your equation A = PDP^(-1) will fail when you check it. P has to be invertible, which means the eigenvectors you collected must be linearly independent. If you found enough eigenvectors for each eigenvalue, they are automatically independent across different eigenvalues. Within the same eigenvalue, you need to verify independence yourself, though in practice the row reduction process usually makes that obvious.
Step Four: Verify By Multiplication
Multiply P * D * P^(-1) and confirm you get A back. Or equivalently, check that A * P = P * D. The second form is computationally cheaper because you do not need to invert P. Just multiply A by each column of P and confirm you get lambda times that column. This catches errors in P construction faster than anything else. Diagonalization turns matrix powers into trivial computations. A^n = P * D^n * P^(-1), and D^n is just each diagonal entry raised to the nth power. This is essential for solving systems of linear differential equations, computing matrix functions, and analyzing discrete dynamical systems. A recurrence relation like x(k+1) = A * x(k) becomes decoupled in the eigenbasis, which is why population models and Markov chain steady states all come back to this. In practice, I use diagonalization when I need A^n for large n or when I am computing e^(At). Without it, you are doing repeated matrix multiplication, which scales poorly. With it, the hard part is the eigenvector computation and the one matrix inversion, and everything after that is cheap.
When Diagonalization Fails And What To Do Instead
Not every matrix is diagonalizable. Defective matrices, which have geometric multiplicity strictly less than algebraic multiplicity for at least one eigenvalue, cannot be diagonalized over the reals or complex numbers. Common examples include matrices with a nonzero nilpotent part, like a Jordan block [[2, 1], [0, 2]]. When diagonalization fails, the Jordan normal form is the fallback, but it is numerically unstable for floating-point computation. If you are working with real data or doing actual engineering work, singular value decomposition or Schur decomposition is usually more appropriate. SVD gives you A = U * Sigma * V^T, which works for any matrix and is the backbone of least squares, low-rank approximation, and regularization. Schur decomposition gives you A = Q * T * Q^T where T is upper triangular, which is close enough to diagonal for most numerical purposes and is what most numerical libraries actually compute internally. Another failure mode is when the matrix is symmetric but you are working in a context that requires orthogonal diagonalization. In that case, you need to normalize your eigenvectors so that P becomes an orthogonal matrix, which means P^(-1) = P^T. This simplifies the computation significantly and is guaranteed to work for any real symmetric matrix by the spectral theorem. I always double-check symmetry before proceeding because it saves me from having to invert a general matrix later.

Practical Notes From Experience
Symbolic computation tools like SymPy or Mathematica can handle the eigenvalue and eigenvector computation exactly, which avoids roundoff issues. But they can be slow on larger matrices and sometimes return eigenvectors in weird forms that require manual simplification. For routine work, I compute eigenvalues numerically with NumPy or MATLAB, then verify the eigenvector equations by substitution to catch any conditioning problems. Conditioning is a real concern. If your eigenvalues are close together or your eigenvectors are nearly linearly dependent, small perturbations in the matrix entries can cause large changes in the eigenvectors. This is measured by the eigenvector condition number, and when it is large, the diagonalization is numerically unreliable. In those cases, stick with Schur or SVD. Also worth noting: diagonalization over the reals requires all eigenvalues to be real. If your characteristic polynomial has complex roots, you can still diagonalize over the complex numbers, but the resulting P and D will have complex entries. For real-valued applications, this sometimes means keeping complex conjugate pairs together and using a real block diagonal form instead, though that is a separate topic.
The entire process from start to finish on a typical 3x3 homework problem takes about ten to fifteen minutes if you know what you are doing and twelve to twenty minutes if you are being careful. A 4x4 with distinct eigenvalues might take twenty to forty minutes by hand, which is why nobody does those by hand anymore. The structure is the same regardless of size, but the bookkeeping gets ugly fast.