generalized Fourier series
The ordinary Fourier series writes a function as a sum of sines and cosines. But sines and cosines are just one orthogonal basis among many — the one that fits a string fixed at both ends. A generalized Fourier series is the same idea with the basis swapped out: expand a function as an infinite sum over the eigenfunctions of whatever Sturm-Liouville problem your geometry produced. Cylinders ask for Bessel functions, spheres for Legendre polynomials, an oscillator for Hermite functions; the bookkeeping is identical.
Write f(x) = sum over n of c_n y_n(x), where {y_n} are the orthogonal eigenfunctions and {c_n} are the generalized Fourier coefficients. Orthogonality with the weight w hands you each coefficient by a single projection: c_n = (integral from a to b of f(x) y_n(x) w(x) dx) divided by (integral from a to b of y_n(x)^2 w(x) dx). The numerator measures how much of mode n is present in f; the denominator (the squared norm) normalizes. If the eigenfunctions are normalized to unit norm, the denominator is 1 and c_n is just the projection. Completeness of the basis is what guarantees the series actually reconstructs f, in the mean-square sense.
This is the practical engine of separation of variables. After separating a PDE you are handed an eigenfunction basis and an initial or boundary profile; expanding that profile in a generalized Fourier series gives the coefficients, and then each mode evolves on its own. The convergence is the honest part to watch: like the classical series, a generalized Fourier series of a discontinuous f converges in mean square but overshoots at jumps, and its termwise derivative may converge much more slowly than the series itself.
To expand f(x) = 1 in Legendre polynomials on [-1, 1], project: c_n = ((2n+1)/2) integral from -1 to 1 of 1 times P_n dx. Only c_0 survives, giving simply f = P_0 = 1 — as it must, since the constant is itself the zeroth Legendre polynomial.
The coefficient formula is one projection per mode; the squared-norm denominator for P_n on [-1, 1] is 2/(2n+1).
The coefficients are unique and easy only because the basis is orthogonal; if you expand in non-orthogonal functions you cannot project mode-by-mode and must instead solve a coupled linear system for all coefficients at once.