Why these equations, and not a thousand others
Across this rung you learned a general craft: feed a linear second-order equation a series, grind out a recurrence relation, and let the coefficients tell you the solution. That craft works on infinitely many equations, most of which have no name and no fame. So why does every textbook stop and name a few — Bessel, Legendre, Hermite, Airy — and lavish whole chapters on them? Because these particular equations are not random: they are exactly what falls out when you try to solve the wave equation, the heat equation, or Laplace's equation in round or spherical regions.
Here is the mechanism. When you solve a partial differential equation by separation of variables, the multi-variable problem splinters into several single-variable ODEs — one per coordinate. In Cartesian coordinates those ODEs are the familiar y'' + lambda y = 0, whose answers are sines and cosines. But change to polar or spherical coordinates and the coordinate factors 1/r and 1/sin(theta) creep into the equation, turning the constant coefficients into variable ones with a singular point sitting at the origin or the pole. The radial part becomes Bessel's equation; the polar-angle part becomes Legendre's. Sines and cosines were the special functions of flat geometry; Bessel and Legendre are the special functions of round geometry.
Legendre: when the series politely stops
Start with the friendlier of the two. The Legendre equation (1 - x^2) y'' - 2x y' + n(n+1) y = 0 has an ordinary point at x = 0, so you can attack it with the straight power-series method from the start of this rung — no Frobenius needed there. Substitute y = sum a_k x^k, collect powers, and the two-term recurrence reads a_(k+2) = a_k · [k(k+1) - n(n+1)] / [(k+1)(k+2)]. Notice the numerator: when k reaches n, the bracket k(k+1) - n(n+1) becomes exactly zero.
That single zero changes everything. Once a_(k+2) = 0, every later coefficient in that chain is zero too, and the infinite series guillotines down to a polynomial of degree n. For each non-negative integer n you get a finite, well-behaved answer — these are the Legendre polynomials P_n(x): P_0 = 1, P_1 = x, P_2 = (3x^2 - 1)/2, and so on. The other independent solution (the chain that does not terminate) blows up at x = +/- 1, the very poles of the sphere, so physics quietly discards it and keeps the polynomial. The requirement "the solution must stay finite at the poles" is what forces n to be a whole number in the first place.
The Legendre polynomials have a gift that single sines and cosines do not flaunt: they are mutually orthogonal on the interval [-1, 1]. The integral of P_m(x) P_n(x) over that interval is zero whenever m is not equal to n. This orthogonality (here with weight 1) is the engine behind expanding an arbitrary function on a sphere into a sum of Legendre pieces — exactly as a Fourier series expands a periodic function into sines and cosines. Each P_n catches one ingredient of the function and the others integrate away to nothing.
Bessel: a genuine Frobenius case at last
Bessel's equation x^2 y'' + x y' + (x^2 - p^2) y = 0 has a regular singular point at x = 0 — a plain power series fails there, and you must reach for the full method of Frobenius you built two guides ago. The indicial equation comes out as r^2 - p^2 = 0, giving roots r = +p and r = -p. This is precisely the showcase the indicial equation was built for: the two roots either differ by a non-integer (the easy case, two clean Frobenius series) or by an integer (the delicate case, where the second solution may drag in a logarithm).
Take the larger root r = +p and turn the Frobenius crank. The recurrence forces all odd coefficients to vanish and ties each even one to the previous, and out comes the Bessel function of the first kind, J_p(x). Written out, J_p(x) = sum over m of [ (-1)^m / (m! · Gamma(m + p + 1)) ] · (x/2)^(2m + p). The factorials in the denominator grow ferociously fast, so the radius of convergence is infinite — J_p is an honest entire-style function, oscillating like a sine but with an amplitude that slowly decays as 1/sqrt(x) and with zeros that are not evenly spaced.
What about the second solution? When p is not an integer, the smaller root r = -p delivers a perfectly good independent series J_(-p)(x), and the general solution is just a combination of the two. But when p is an integer the two roots differ by an integer and the easy second series collapses into a multiple of the first — exactly the snag the indicial-equation guide warned you about. Frobenius then forces a logarithm into the second solution, called the Bessel function of the second kind Y_p(x), which is unbounded at x = 0. The full general solution is y = c1 J_p(x) + c2 Y_p(x), and in a physical disk where the center must stay finite, you keep J_p and throw Y_p away.
A small zoo, and what they share
Bessel and Legendre have cousins, each the named solution of its own variable-coefficient equation, each pulled from physics. Here is the family at a glance.
function equation / setting shows up in -------------------------------------------------------------------- Legendre P_n (1-x^2)y'' - 2x y' + n(n+1)y=0 sphere, gravity/electrostatic potential Bessel J_p x^2 y'' + x y' + (x^2-p^2)y = 0 drumhead, cylinder, wave in a pipe Hermite H_n y'' - 2x y' + 2n y = 0 quantum harmonic oscillator Laguerre L_n x y'' + (1-x)y' + n y = 0 hydrogen atom radial part Chebyshev T_n (1-x^2)y'' - x y' + n^2 y = 0 approximation, numerical methods Airy Ai, Bi y'' - x y = 0 turning points, optics caustics
The pattern is striking once you see it. Hermite, Laguerre and Chebyshev all share Legendre's good fortune: for integer n their series terminate into polynomials, and they are orthogonal — each with its own weight function, e^(-x^2) for Hermite, e^(-x) for Laguerre, 1/sqrt(1-x^2) for Chebyshev. The Airy functions, by contrast, never terminate; they are full power series that behave like decaying waves on one side of x = 0 and growing exponentials on the other, which is why they govern the transition at a turning point. And many of these are special cases of one grand parent, the hypergeometric equation, whose series contains Legendre, Chebyshev and more as particular settings of its parameters.
Orthogonality is the deeper reason they matter
It is tempting to think the point of these functions is to write down closed forms. The deeper payoff is the orthogonality you glimpsed with Legendre, and it is no accident. Each of these equations can be massaged into a common shape — a Sturm-Liouville problem — and a theorem guarantees that the solutions belonging to different eigenvalues are orthogonal with respect to a specific weight. That orthogonality is what lets you expand a general function as a series of special functions and solve for the coefficients one integral at a time, the same divide-and-conquer logic behind Fourier series.
So the chain of ideas closes neatly. A PDE in round coordinates separates into named ODEs; Frobenius and the power-series method solve those ODEs and produce Bessel, Legendre and their cousins; the physical demand of finiteness picks out the well-behaved branch and quantizes the order n; and orthogonality then lets you reassemble the full PDE solution as a series of these pieces, fitting the boundary and initial data coefficient by coefficient. Every tool from this rung is doing real work in that chain.
One honest caveat to carry forward. Naming a function does not mean you can evaluate it by hand: J_2(3.7) or P_5(0.41) come from tables, recurrences, or a library call, never from a tidy formula in elementary functions. And these named functions, glorious as they are, cover only the small set of equations that happen to separate the standard geometries — the vast majority of variable-coefficient ODEs have no special-function answer at all and must be handed to the numerical and qualitative methods of the next rungs. The special functions are a beautiful, finite island of solvability in a much larger sea.