eigenfunction expansion method
Here is the master strategy that ties this whole field to the solving of real boundary-value problems: to solve a differential equation, expand everything in the eigenfunctions of the problem's own operator. Because those eigenfunctions diagonalize the operator — acting on each one just multiplies it by its eigenvalue — a complicated differential equation collapses into a list of trivial, decoupled equations for the coefficients, one per mode.
Concretely, take a boundary-value or initial-value problem governed by a Sturm-Liouville operator L. Find its eigenfunctions y_n and eigenvalues lambda_n. Expand the unknown solution, the source term, and the initial data all in the same basis. Because L y_n = -lambda_n w y_n, the operator becomes multiplication by -lambda_n inside the sum, so each coefficient obeys its own simple equation. For a forced steady problem L u = f, you get coefficient c_n = (projection of f onto y_n) divided by (-lambda_n times the squared norm) — provided no lambda_n is zero (if one is, the Fredholm alternative warns you a solution exists only if f is orthogonal to that mode). For a time-dependent problem like the heat equation, each mode's amplitude decays as e^{-lambda_n t}, so the solution is sum of c_n e^{-lambda_n t} y_n(x).
This is exactly what separation of variables is doing under the hood, and it is also how a Green's function gets built as a sum over eigenfunctions. The method's reach and its limits are both honest: it is powerful and exact for linear problems on domains where you can actually find the eigenfunctions (intervals, disks, spheres, boxes in nice coordinates), but it does not help with nonlinear equations or irregular domains where no clean eigenbasis exists, and like every such expansion it inherits mean-square convergence and the Gibbs overshoot at discontinuities.
Solve the heat equation u_t = u_xx on [0, L] with u = 0 at both ends and initial profile u(x, 0) = f(x). Expand f in the eigenfunctions sin(n pi x / L): then u(x, t) = sum of b_n e^{-(n pi/L)^2 t} sin(n pi x / L), each mode decaying at its own rate, the lowest surviving longest.
Each eigenmode evolves independently; high modes (large lambda_n) vanish fastest, which is why diffusion smooths sharp features first.
The method is exact and clean only when a complete eigenbasis is actually available; for irregular domains or nonlinear equations there may be no usable eigenfunctions, and a zero eigenvalue makes a steady forced problem solvable only under the Fredholm-alternative orthogonality condition.