From dividing the function to inverting the operator
Look back at the path you have climbed in this rung. You learned to pose a boundary value problem, you found that the homogeneous version has special values lambda where nontrivial solutions exist — the eigenvalues and their eigenfunctions — and you used the magic of orthogonality to write any function as an eigenfunction expansion. That expansion is a beautiful solver: to solve L y = f, you expand both y and f in eigenfunctions and match coefficient by coefficient. But notice what it costs you. You must first find every eigenfunction, then compute an infinite stack of integrals, then resum the series. It is divide-and-conquer taken to the level of the unknown function itself.
A Green's function asks a sharper question. Instead of dividing up the answer y, it divides up the differential operator's *job*. The boundary value problem reads L y = f, where L is a linear differential operator like L = -d^2/dx^2 plus boundary conditions glued on. Now L is linear — it obeys superposition — so it behaves like a matrix acting on functions. And what do you do to a matrix equation A x = f? You invert: x = A^(-1) f. The whole game of this guide is to build the inverse of L, an object that takes any f and hands you the solution y directly, in one integral, no eigenfunction hunt required.
The unit poke: response to a delta
How do you invert a matrix conceptually? You see how it acts on each basis vector e_j and assemble the columns of A^(-1). The continuous version of a basis vector — a spike concentrated entirely at one point — is the Dirac delta delta(x - s): zero everywhere except at x = s, yet with total integral 1. Think of it as a unit hammer-tap delivered at the single location s. The Green's function G(x, s) is, quite literally, the system's response measured at x when you tap it at s. It solves L G = delta(x - s), together with the same boundary conditions the original problem carried.
Here is the payoff, and it is the whole idea in one line. Any forcing f(x) is just a continuous pile of tiny pokes: f(x) = integral of f(s) delta(x - s) ds. By superposition, the response to the pile is the pile of responses. So the solution of L y = f is simply y(x) = integral of G(x, s) f(s) ds. The forcing f is the recipe of how hard to tap at each location s; G(x, s) is the response to a single tap; and the integral adds up all those responses, weighted by f, to reconstruct the full answer. The Green's function carries all the difficulty — once you have it, every forcing f is solved by the same one integral.
Building G by hand: a continuous, kinked thing
What does G(x, s) actually look like? Away from the tap point x = s the delta is zero, so G satisfies the *homogeneous* equation L G = 0 there. So on the left of s, G is some homogeneous solution chosen to obey the left boundary condition; on the right of s, another homogeneous solution chosen to obey the right boundary condition. The two pieces meet at x = s, and the delta lives entirely in how they are joined. This is why a Green's function always comes in two halves — a function defined piecewise, with the seam at x = s.
Two matching conditions stitch the halves at the seam, and they are the soul of the construction. First, G must be continuous across x = s — no jump in the value itself. Second, its slope must have a precise jump: the derivative dG/dx leaps by a fixed amount as you cross s. Why a kink and not a smooth join? Because the second derivative of a kink is a spike — exactly the delta we are trying to produce. A pure corner in G is the only way for -d^2G/dx^2 to deliver a delta of unit strength. For the operator -d^2/dx^2 the jump in slope is exactly -1; for a general Sturm-Liouville operator it is set by the leading coefficient p(s).
Build G(x,s) for -y'' = f(x), y(0)=0, y(1)=0 on [0,1]:
homogeneous solutions of -y'' = 0 are 1 and x
x < s : G = a*x (kills y(0)=0 since G(0)=0)
x > s : G = b*(1 - x) (kills y(1)=0 since G(1)=0)
match at x = s:
continuity : a*s = b*(1 - s)
slope jump : b*(-1) - a = -1 ( G' jumps by -1 )
solve -> a = 1 - s , b = s
G(x,s) = x(1 - s) for x <= s
s(1 - x) for x >= s ( note: G(x,s) = G(s,x) )Symmetry, self-adjointness, and the spectral formula
Stare at the answer above: G(x, s) = G(s, x). The response felt at x to a tap at s equals the response felt at s to a tap at x. This symmetry is not a lucky accident of the example — it is the continuous echo of the fact that a symmetric matrix has a symmetric inverse. It holds precisely when L is written in self-adjoint form, the very form a Sturm-Liouville problem is engineered to wear. The same self-adjointness that earlier guaranteed real eigenvalues and orthogonal eigenfunctions now guarantees a symmetric Green's function. The threads of this rung are all the same thread.
And here the two big solvers of this rung finally shake hands. Take the Green's function and expand *it* in the eigenfunctions phi_n of L. Because L phi_n = lambda_n phi_n, inverting L just divides each piece by its eigenvalue, and you get the spectral formula: G(x, s) = sum over n of phi_n(x) phi_n(s) / lambda_n. Read it slowly — it says the inverse operator is built by taking each eigenfunction, dividing by its eigenvalue, and summing. That is the function-world copy of A^(-1) = sum of (eigenvector)(eigenvector)^T / lambda for a symmetric matrix. The eigenfunction expansion and the Green's function were never rivals; they are two faces of inverting the same operator.
Why bother — and the honest fine print
What have you actually bought? A Green's function is computed *once* for a given operator and boundary conditions, and then it solves L y = f for every forcing f by the single integral y(x) = integral of G(x, s) f(s) ds. Change the load on your beam, the heat source along your rod, the charge in your region — re-run one integral, never re-solve the differential equation. This is exactly why physics and engineering reach for Green's functions: they are the reusable inverse of a system, the response to a unit source that, once known, answers every question about that system. They also slot neatly beside variation of parameters, which builds essentially this same kernel from a fundamental pair of solutions.
Now the honest limits, because elegance is not omnipotence. The clean construction above leans on linearity: a Green's function exists for linear operators only, and the moment your equation is nonlinear, superposition collapses and the whole edifice — delta decomposition, summing of responses, the inverse operator — falls with it. There is no Green's function for a nonlinear boundary value problem; those need iteration or numerics. The method also assumes you can actually solve the homogeneous equation to build the two pieces of G, which for variable-coefficient operators may itself have no closed form.
Two more cautions worth carrying. First, the delta function is not a function in the ordinary sense — it is a distribution, a rigorous object only once you treat it as a recipe inside an integral; the hand-wavy "infinitely tall spike" picture is a useful lie that the integral definition makes honest. Second, the Green's function depends on the boundary conditions as much as on the operator: change y(1) = 0 to y'(1) = 0 and you must rebuild G from scratch — it is the inverse of *the whole problem*, side conditions and all, not of the bare operator alone. Keep both in mind and the Green's function becomes the most powerful single tool in your boundary-value kit.