Formula & theorem reference

This result is proved by constructing a sequence of successive approximations y_0(x), y_1(x), y_2(x), \dots (Picard iterates), each defined from the last by y_{n+1}(x) = y_0 + \int_{x_0}^{x} f(t, y_n(t))\,dt, and showing this sequence converges to a genuine solution — the continuity of f alone guarantees a solution exists (this weaker half is Peano's theorem). The extra hypothesis, continuity of \partial f/\partial y, is exactly what makes the map y_n \mapsto y_{n+1} a contraction on a small enough interval: two different starting guesses get pulled closer together with each iteration, at a rate bounded by a Lipschitz constant derived from \partial f/\partial y. A genuine second solution through the same point would be a second fixed point of that contraction, which a contraction mapping can never have. Neither half of this argument produces a formula — both are pure existence arguments, which is exactly why this theorem answers "is there a solution, and only one" without ever constructing it explicitly.

If the equation is exact, some potential function F(x,y) satisfies F_x = M and F_y = N. Since M and N have continuous partials, Clairaut's theorem on mixed partials gives M_y = F_{xy} = F_{yx} = N_x, proving the condition is necessary. For sufficiency, construct F directly: integrate M with respect to x, F(x,y) = \int M(x,y)\,dx + h(y), where h(y) is an unknown function of y alone. Differentiating this candidate with respect to y and setting the result equal to N isolates h'(y); the hypothesis M_y = N_x is exactly what guarantees h'(y) comes out depending on y alone (never x), so it can be integrated to find h(y) and complete F.

Multiplying by \mu(x) turns the equation into \mu M\,dx + \mu N\,dy = 0, exact exactly when (\mu M)_y = (\mu N)_x. Since \mu depends only on x, (\mu M)_y = \mu M_y, while (\mu N)_x = \mu' N + \mu N_x by the product rule. Setting these equal and dividing by \mu N gives \dfrac{\mu'}{\mu} = \dfrac{M_y - N_x}{N} — a separable equation for \mu exactly when the right side depends on x alone, and it integrates directly to \ln\mu = \int\frac{M_y-N_x}{N}\,dx, giving the stated formula. The \mu(y) case follows by the same argument with the roles of M_y and N_x exchanged.

Divide the original equation by y^n: y^{-n}\dfrac{dy}{dx} + P(x)y^{1-n} = Q(x). With v = y^{1-n}, the chain rule gives \dfrac{dv}{dx} = (1-n)y^{-n}\dfrac{dy}{dx}, so y^{-n}\dfrac{dy}{dx} = \dfrac{1}{1-n}\dfrac{dv}{dx}. Substituting both this and y^{1-n}=v into the divided equation gives \dfrac{1}{1-n}\dfrac{dv}{dx} + P(x)v = Q(x), and multiplying through by 1-n produces exactly the stated linear equation in v.

With y = vx, the product rule gives \dfrac{dy}{dx} = v + x\dfrac{dv}{dx}. Substituting into \dfrac{dy}{dx} = g(y/x) = g(v) gives v + x\dfrac{dv}{dx} = g(v), and isolating the derivative term produces x\dfrac{dv}{dx} = g(v) - v — separable in v and x, since the right side depends only on v.

Near the critical parameter value, a smooth f(y,c) with a double root at (y^*, c^*) (so f=0 and \partial f/\partial y=0 there, but \partial f/\partial c \ne 0 and \partial^2 f/\partial y^2 \ne 0) has a Taylor expansion in (y-y^*) and (c-c^*) whose leading surviving terms are exactly linear in (c-c^*) and quadratic in (y-y^*) — the same shape as c-y^2 after an affine change of variables and a rescaling of c. This is why the equilibrium count and stability pattern (two, one, zero) is universal for this kind of collision: it depends only on which derivatives vanish, not on the specific formula f started as.

A method built by matching a Taylor expansion through the h^k term has local truncation error O(h^{k+1}) — one power of h higher than what it matched, since the first term it fails to match is the source of the error. Reaching a fixed endpoint from a fixed starting point takes n = (\text{interval length})/h steps, so the global error is (local error per step) \times (number of steps) \sim h^{k+1}\cdot(1/h) = h^k — one power of h lower than the local error. This is why Euler's method (matches through h^1, local error O(h^2)) is global order 1, while Heun's method (matches through h^2, local error O(h^3)) is global order 2, and RK4 (matches through h^4, local error O(h^5)) is global order 4: reading a method's global order directly off its local error requires subtracting exactly one power of h.

k_1 is an ordinary Euler step's slope, used to produce a rough predicted point \tilde y at x_{n+1}. k_2 is the slope f actually has there. Averaging the two slopes, (k_1+k_2)/2, approximates the slope a trapezoid connecting (x_n,y_n) to (x_{n+1}, \tilde y) would use — the same idea the trapezoid rule uses to approximate a definite integral, applied one step at a time to a differential equation instead of a fixed integrand. Because this averaged slope matches the Taylor expansion of the true solution through the h^2 term (Euler's single slope only matches it through the h^1 term), the local truncation error is O(h^3) per step, and the accumulated global error over a fixed interval is O(h^2) — one full order better than Euler's method.

k_1 is the slope at the left endpoint; k_2 and k_3 are two independent estimates of the slope at the midpoint, the second refining the first by using k_2's own predicted y-value instead of k_1's; k_4 is the slope at the right endpoint, predicted using the midpoint slope k_3. Weighting these 1:2:2:1 and dividing by 6 reproduces Simpson's rule's own weighting for approximating an integral by samples at both endpoints and (doubly) the midpoint. Carrying this weighted average through a Taylor expansion of the true solution shows it matches every term through h^4, giving a local truncation error of O(h^5) per step and a global error of O(h^4) — two full orders better than Euler's method, at the cost of four function evaluations per step instead of one.

Suppose c_1y_1 + c_2y_2 = 0 for all x, not both c_1, c_2 zero (the definition of dependence). Differentiating gives c_1y_1' + c_2y_2' = 0 too, so (c_1, c_2) solves the linear system \begin{pmatrix} y_1 & y_2 \\ y_1' & y_2' \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} = 0 nontrivially at every x — which forces that matrix's determinant, W(x), to vanish identically. Contrapositively, W nonzero anywhere rules out dependence, proving independence. The "all-or-nothing" half (one zero forces all zero) follows from Abel's formula: W(x) = W(x_0)e^{-\int_{x_0}^{x} p(t)\,dt}, and an exponential is never zero, so W(x_0)=0 forces W(x)=0 for every x.

Differentiate W = y_1y_2' - y_1'y_2 directly: W' = y_1y_2'' - y_1''y_2 (the y_1'y_2' terms cancel). Since both y_i solve the equation, y_i'' = -p y_i' - q y_i. Substituting: W' = y_1(-py_2'-qy_2) - (-py_1'-qy_1)y_2 = -p(y_1y_2'-y_1'y_2) = -pW. So W itself solves the first-order linear equation W' = -pW, whose solution is exactly W(x) = W(x_0)e^{-\int_{x_0}^x p(t)\,dt} — the ordinary exponential-decay formula from Chapter 2, applied to W rather than to y_1 or y_2 themselves.

Let y = c_1y_1 + c_2y_2. Then y' = c_1y_1' + c_2y_2' and y'' = c_1y_1'' + c_2y_2'', since differentiation is linear. Substituting into the equation: y'' + py' + qy = c_1(y_1''+py_1'+qy_1) + c_2(y_2''+py_2'+qy_2) = c_1(0) + c_2(0) = 0, using that y_1, y_2 each already solve the equation. The result holds for a homogeneous equation specifically because its right side is 0 — c_1(0)+c_2(0) collapses to 0 only when both individual right sides are 0, which is exactly why superposition fails for a nonhomogeneous equation's own solutions (later chapters).

Substitute x = i\theta into the Maclaurin series for e^x and use i^2=-1, i^3=-i, i^4=1 (the pattern repeats every four powers): e^{i\theta} = \sum_{n=0}^{\infty}\frac{(i\theta)^n}{n!} = 1 + i\theta - \frac{\theta^2}{2!} - i\frac{\theta^3}{3!} + \frac{\theta^4}{4!} + i\frac{\theta^5}{5!} - \cdots Collecting the terms with no factor of i and the terms with one factor of i into two separate sums gives e^{i\theta} = \underbrace{\left(1 - \frac{\theta^2}{2!} + \frac{\theta^4}{4!} - \cdots\right)}_{\cos\theta} + i\underbrace{\left(\theta - \frac{\theta^3}{3!} + \frac{\theta^5}{5!} - \cdots\right)}_{\sin\theta}. The first parenthesis is exactly the Maclaurin series for \cos\theta and the second is exactly the Maclaurin series for \sin\theta, so e^{i\theta} = \cos\theta + i\sin\theta.

Substituting y=x^r gives y'=rx^{r-1}, y''=r(r-1)x^{r-2}, so ax^2y''+bxy'+cy = \big[ar(r-1)+br+c\big]x^r — this vanishes exactly when r solves ar(r-1)+br+c=0, the equation's own characteristic equation, structurally identical to the constant-coefficient case with x^r standing in for e^{rx}. The substitution t=\ln x turns the Cauchy-Euler equation into an ordinary constant-coefficient equation in t directly (since x\frac{dy}{dx}=\frac{dy}{dt}), which is why \ln x appears in the repeated-root and complex cases exactly where x (a bare t) appeared in the constant-coefficient formulas.

This is Chapter 7's constant-coefficient classification applied verbatim to m, c, k in place of a, b, c — the discriminant's sign decides distinct real, repeated, or complex-conjugate roots, and c\ge0 together with m,k>0 rules out any root with positive real part, so every case decays or oscillates in place rather than growing. Dividing c^2-4mk<0 through by 4m^2 gives (c/(2m))^2<k/m=\omega_0^2, i.e. \zeta^2<1 — the discriminant test and the \zeta-versus-1 test are the same inequality in different units.

Kirchhoff's voltage law around a series RLC loop sums the voltage drops — L\,dI/dt across the inductor, RI across the resistor, q/C across the capacitor — equal to the source voltage E(t). Substituting I=q' gives Lq''+Rq'+q/C=E(t) directly, term for term the same linear differential equation as mx''+cx'+kx=F(t), so every theorem already proved for one applies to the other without re-derivation.

An undetermined-coefficients trial works by substituting y_p into ay''+by'+cy and matching coefficients against g(x) — this only produces a solvable system when the trial's own functional form is closed under differentiation (each derivative of a polynomial is a lower-degree polynomial, each derivative of e^{\alpha x} is a multiple of itself, each derivative of \cos/\sin(\beta x) stays inside the same two-dimensional span). If the naive trial already solves the homogeneous equation, the linear operator ay''+by'+cy annihilates it entirely — substituting it produces 0=g(x), an unsolvable equation regardless of the coefficients chosen. Multiplying by x^s breaks the trial out of the homogeneous solution space (the same fix reduction of order uses for a repeated root) while staying closed under differentiation, restoring a solvable system.

Seek y_p=u_1y_1+u_2y_2 where u_1,u_2 are functions (not constants) to be determined — hence "variation" of the parameters C_1,C_2 in the homogeneous general solution. Differentiating and imposing the convenience condition u_1'y_1+u_2'y_2=0 (which removes second derivatives of u_1,u_2 from the result) gives y_p'=u_1y_1'+u_2y_2', and then y_p''=u_1'y_1'+u_2'y_2'+u_1y_1''+u_2y_2''. Substituting into the original equation and using that y_1,y_2 each already solve the homogeneous equation collapses everything except u_1'y_1'+u_2'y_2'=g. Paired with the convenience condition u_1'y_1+u_2'y_2=0, this is a linear system in u_1',u_2' with coefficient determinant W, solved by Cramer's rule as u_1'=-y_2g/W and u_2'=y_1g/W. Integrating and substituting back into y_p=u_1y_1+u_2y_2 produces the boxed formula. Because W\ne0 for any fundamental set (Chapter 6's Wronskian test), the formula applies to every constant-coefficient or variable-coefficient linear equation for which a fundamental set is already known — unlike undetermined coefficients, it places no restriction on g(x)'s form.

Let y = c_1y_1 + c_2y_2. Then y' = c_1y_1' + c_2y_2' and y'' = c_1y_1'' + c_2y_2'', since differentiation is linear. Substituting into the equation: y'' + py' + qy = c_1(y_1''+py_1'+qy_1) + c_2(y_2''+py_2'+qy_2) = c_1(0) + c_2(0) = 0, using that y_1, y_2 each already solve the equation. The result holds for a homogeneous equation specifically because its right side is 0 — c_1(0)+c_2(0) collapses to 0 only when both individual right sides are 0, which is exactly why superposition fails for a nonhomogeneous equation's own solutions (later chapters).

Seek y_p=u_1y_1+u_2y_2 where u_1,u_2 are functions (not constants) to be determined — hence "variation" of the parameters C_1,C_2 in the homogeneous general solution. Differentiating and imposing the convenience condition u_1'y_1+u_2'y_2=0 (which removes second derivatives of u_1,u_2 from the result) gives y_p'=u_1y_1'+u_2y_2', and then y_p''=u_1'y_1'+u_2'y_2'+u_1y_1''+u_2y_2''. Substituting into the original equation and using that y_1,y_2 each already solve the homogeneous equation collapses everything except u_1'y_1'+u_2'y_2'=g. Paired with the convenience condition u_1'y_1+u_2'y_2=0, this is a linear system in u_1',u_2' with coefficient determinant W, solved by Cramer's rule as u_1'=-y_2g/W and u_2'=y_1g/W. Integrating and substituting back into y_p=u_1y_1+u_2y_2 produces the boxed formula. Because W\ne0 for any fundamental set (Chapter 6's Wronskian test), the formula applies to every constant-coefficient or variable-coefficient linear equation for which a fundamental set is already known — unlike undetermined coefficients, it places no restriction on g(x)'s form.

Integrate the defining transform of f' by parts, with u = e^{-st}, dv = f'(t)\,dt: \mathcal{L}\{f'(t)\} = \int_0^\infty f'(t)e^{-st}\,dt = \left[f(t)e^{-st}\right]_0^\infty + s\int_0^\infty f(t)e^{-st}\,dt. The boundary term at t=0 contributes -f(0); the term at t\to\infty vanishes for every s in the transform's region of convergence, since e^{-st} decays faster than f can grow there. What remains is sF(s) - f(0).

For f'', apply this same rule to g = f': \mathcal{L}\{g'(t)\} = s\mathcal{L}\{g(t)\} - g(0) = s\big(sF(s)-f(0)\big) - f'(0) = s^2F(s) - sf(0) - f'(0). Each further derivative repeats the same integration-by-parts step once more, producing one more factor of s and one more initial-condition term — the pattern the general nth- order statement packages. This is exactly what makes the transform solve an IVP directly (§12): every derivative in the equation becomes a power of s times F(s), with the initial conditions built in from the start rather than applied afterward.

Substitute directly into the defining integral: \mathcal{L}\{e^{at}f(t)\} = \int_0^\infty e^{at}f(t)e^{-st}\,dt = \int_0^\infty f(t)e^{-(s-a)t}\,dt = F(s-a), recognizing the last integral as F evaluated at s-a instead of s. The shifted region of convergence follows because the original integral needed s - a > \alpha for convergence, i.e. s > \alpha + a.

Split the defining integral at t=c, where u_c(t)f(t-c) is 0 for t<c and f(t-c) for t \ge c: \mathcal{L}\{u_c(t)f(t-c)\} = \int_c^\infty f(t-c)e^{-st}\,dt. Substitute \tau = t - c, so t = \tau + c and the lower limit becomes \tau = 0: \int_0^\infty f(\tau)e^{-s(\tau+c)}\,d\tau = e^{-sc}\int_0^\infty f(\tau)e^{-s\tau}\,d\tau = e^{-cs}F(s). This is the algebraic reason a piecewise-defined forcing function built from unit steps (§9) transforms termwise: each delayed piece contributes its own e^{-cs} factor to an otherwise ordinary transform.

Write out the transform of the convolution as a double integral and swap the order of integration: \mathcal{L}\{(f*g)(t)\} = \int_0^\infty\!\!\int_0^t f(\tau)g(t-\tau)\,d\tau\,e^{-st}\,dt = \int_0^\infty\!\!\int_\tau^\infty f(\tau)g(t-\tau)e^{-st}\,dt\,d\tau, where the region \{0 \le \tau \le t\} is now described by t ranging over [\tau, \infty) for each fixed \tau \ge 0. Substitute u = t - \tau in the inner integral (t = u+\tau, dt = du, lower limit u=0): \int_0^\infty f(\tau)\left(\int_0^\infty g(u)e^{-s(u+\tau)}\,du\right)d\tau = \int_0^\infty f(\tau)e^{-s\tau}\,d\tau \int_0^\infty g(u)e^{-su}\,du = F(s)G(s), since the inner integral factors as e^{-s\tau} times G(s), leaving exactly F(s) times G(s) once the outer integral is recognized.

Substitute \mathbf{x}(t)=\mathbf{v}e^{\lambda t} into \mathbf{x}'=A\mathbf{x}: the left side is \mathbf{x}'(t)=\lambda\mathbf{v}e^{\lambda t}, and the right side is A\mathbf{x}(t)=A\mathbf{v}e^{\lambda t}. These are equal for every t exactly when \lambda\mathbf{v}=A\mathbf{v} — the defining equation of an eigenpair. So every eigenpair (\lambda,\mathbf{v}) of A produces one genuine solution this way. Each such solution solves a linear system, so by the superposition principle (the systems analogue of superposition-principle) any linear combination of these n solutions is again a solution; when the n eigenvectors are linearly independent, the resulting solutions are linearly independent functions of t (their values at t=0 are exactly \mathbf{v}_1,\ldots, \mathbf{v}_n, already independent), so their combination is the full n-parameter general solution, not just a solution.

Differentiate the series term by term (justified since it converges for every A and t, by the same ratio-test argument that shows e^x's series converges everywhere): \frac{d}{dt}e^{At} = A + A^2t + \frac{A^3t^2}{2!} + \cdots = A\left(I+At+\frac{(At)^2}{2!}+\cdots\right) = Ae^{At}. So \mathbf{x}(t)=e^{At}\mathbf{x}_0 satisfies \mathbf{x}'(t)=Ae^{At}\mathbf{x}_0=A\mathbf{x}(t), and at t=0 the series collapses to e^{A\cdot0}=I, giving \mathbf{x}(0)=\mathbf{x}_0 exactly as required. Uniqueness follows from Picard–Lindelöf (ch.01): A\mathbf{x} is linear, hence Lipschitz in \mathbf{x}, so the IVP has exactly one solution, and this formula is it.

For the fundamental-matrix identity: any fundamental matrix \Phi(t) satisfies \Phi'(t)=A\Phi(t) column by column, the same equation e^{At} satisfies, so \Phi(t)\Phi(0)^{-1} does too (constant matrices factor through differentiation), and at t=0 it equals \Phi(0)\Phi(0)^{-1}=I — the same initial condition e^{At} has. Two solutions of the same matrix IVP that agree at t=0 are equal for all t, by the same uniqueness argument applied column by column.

(\Leftarrow) n independent eigenvectors give A = PDP^{-1}. Suppose v_1, \dots, v_n are linearly independent eigenvectors of A, with Av_i = \lambda_i v_i. Let P be the matrix with these vectors as its columns, and let D = \operatorname{diag}(\lambda_1, \dots, \lambda_n). Multiplying A into P column by column, AP = A\begin{bmatrix} v_1 & \cdots & v_n \end{bmatrix} = \begin{bmatrix} Av_1 & \cdots & Av_n \end{bmatrix} = \begin{bmatrix} \lambda_1 v_1 & \cdots & \lambda_n v_n \end{bmatrix}. Multiplying P by D on the right scales each column of P by the matching diagonal entry of D, producing the very same matrix: PD = \begin{bmatrix} v_1 & \cdots & v_n \end{bmatrix} \operatorname{diag}(\lambda_1,\dots,\lambda_n) = \begin{bmatrix} \lambda_1 v_1 & \cdots & \lambda_n v_n \end{bmatrix}. So AP = PD. Because P's columns are n linearly independent vectors in \mathbb{R}^n, P is invertible, and right-multiplying both sides by P^{-1} gives A = PDP^{-1}.

(\Rightarrow) A = PDP^{-1} gives n independent eigenvectors. Suppose A = PDP^{-1} with P invertible and D = \operatorname{diag}(\lambda_1, \dots, \lambda_n). Then AP = PD, and reading that equation column by column exactly as above shows Av_i = \lambda_i v_i, where v_i is the i-th column of P. Each v_i is nonzero — a zero column would make P singular — so v_1,\dots,v_n are eigenvectors of A. Since P is invertible, its columns are linearly independent, which produces the needed n independent eigenvectors.

Sanity check: take A = \begin{bmatrix} 4 & 1 \\ 0 & 3 \end{bmatrix}, upper triangular, so its eigenvalues sit on the diagonal: \lambda_1 = 4, \lambda_2 = 3. Solving (A - 4I)v = 0 gives eigenvector v_1 = (1,0), and (A-3I)v = 0 gives v_2 = (1,-1). These two are not scalar multiples of each other, so they're independent, and the theorem guarantees A = PDP^{-1} with P = \begin{bmatrix} 1 & 1 \\ 0 & -1 \end{bmatrix} and D = \begin{bmatrix} 4 & 0 \\ 0 & 3 \end{bmatrix} — which a direct multiplication PDP^{-1} confirms reproduces A.

The characteristic equation of A is \lambda^2-T\lambda+D=0, so by the quadratic formula \lambda=\dfrac{T\pm\sqrt{T^2-4D}}{2} — exactly the discriminant T^2-4D that separates real (T^2\ge4D) from complex (T^2<4D) eigenvalues, matching the table's parabola boundary.

If D<0, the product of the eigenvalues \lambda_1\lambda_2=D is negative, so the eigenvalues are real (a negative product under a real square root only happens when T^2-4D>0, since D<0 forces -4D>0) with opposite sign — one stretches trajectories outward, the other pulls them inward, the saddle picture, regardless of T.

If D>0 and T^2>4D (real, same-sign eigenvalues, since their product D is positive), both eigenvalues share the sign of their sum T=\lambda_1+\lambda_2: both negative (stable node) when T<0, both positive (unstable node) when T>0. At T^2=4D the two real eigenvalues coincide — the repeated-eigenvalue, degenerate-node boundary.

If D>0 and T^2<4D, the eigenvalues are a complex conjugate pair \alpha\pm i\beta with \alpha=T/2 (their real part is half their sum) and \beta=\sqrt{4D-T^2}/2. The real part \alpha controls growth/decay of the spiral's envelope e^{\alpha t}: \alpha<0 spirals inward (stable), \alpha>0 spirals outward (unstable), and \alpha=0 — exactly T=0 — leaves a pure rotation of constant radius, the center case, neither growing nor decaying.

Finally D=0 means \det(A)=0, so A has a nontrivial null space — every vector in it is an equilibrium, not an isolated one, the degenerate case the table's last row names.

Write L[y]=(py')'+qy for the Sturm–Liouville operator, so the equation reads L[y]=-\lambda\sigma y. The whole proof rests on one algebraic identity, true for any two twice-differentiable u,v on [a,b]: vL[u]-uL[v] = v(pu')'-u(pv')' = \frac{d}{dx}\Big[p\,(vu'-uv')\Big], which expands by the product rule and the two terms pv'u'-pu'v' cancelling. Integrating both sides from a to b gives Lagrange's identity: \int_a^b\big(vL[u]-uL[v]\big)\,dx = \Big[p\,(vu'-uv')\Big]_a^b. \quad(\ast)

The boundary term vanishes for any u,v satisfying the same separated conditions. At x=a, both (u(a),u'(a)) and (v(a),v'(a)) solve the same linear equation \alpha_1w+\alpha_2w'=0 in the plane — a single line through the origin, since (\alpha_1,\alpha_2)\neq(0,0). Two points on the same line through the origin are proportional (or one is the origin), so the 2\times2 determinant built from them vanishes: v(a)u'(a)-u(a)v'(a)=0. The identical argument at x=b (using \beta_1,\beta_2) kills the term there too, so \big[p(vu'-uv')\big]_a^b=0 whenever u and v both satisfy the problem's boundary conditions.

Reality. Suppose \lambda is an eigenvalue, possibly complex, with eigenfunction y (possibly complex-valued): L[y]=-\lambda\sigma y. Since p,q,\sigma are real-valued, conjugating the whole equation gives L[\bar y]=-\bar\lambda\sigma\bar y. The boundary conditions have real coefficients \alpha_1,\alpha_2,\beta_1,\beta_2, so \bar y satisfies exactly the same conditions y does — the vanishing argument above applies to the pair u=y, v=\bar y. Apply (\ast): \int_a^b\big(\bar yL[y]-yL[\bar y]\big)\,dx = 0 \quad\Longrightarrow\quad \int_a^b\big(-\lambda\,\sigma|y|^2+\bar\lambda\,\sigma|y|^2\big)\,dx = 0 \quad\Longrightarrow\quad (\bar\lambda-\lambda)\int_a^b\sigma|y|^2\,dx = 0. An eigenfunction is not identically zero and \sigma(x)>0 on [a,b], so \int_a^b\sigma|y|^2\,dx>0. Dividing forces \bar\lambda=\lambda: every eigenvalue is real. (Once \lambda is known real, the real and imaginary parts of y each solve the same real equation with the same real \lambda, so a real eigenfunction can always be chosen — used freely from here on.)

Orthogonality. Now take y_m,y_n to be real eigenfunctions for distinct real eigenvalues \lambda_m\neq\lambda_n: L[y_m]=-\lambda_m\sigma y_m and L[y_n]=-\lambda_n\sigma y_n. Apply (\ast) to u=y_m, v=y_n — both satisfy the problem's boundary conditions, so the right-hand side is again 0: \int_a^b\big(y_nL[y_m]-y_mL[y_n]\big)\,dx = 0 \quad\Longrightarrow\quad (\lambda_n-\lambda_m)\int_a^b\sigma\,y_my_n\,dx = 0. Since \lambda_m\neq\lambda_n, dividing gives \int_a^b\sigma(x)y_m(x)y_n(x)\,dx=0 — orthogonality, weighted by \sigma.

Sanity check. Take p\equiv1, q\equiv0, \sigma\equiv1, [a,b]=[0,L], with y(0)=y(L)=0 — this is exactly the eigenvalue problem y''+\lambda y=0 behind the Fourier sine series, whose eigenfunctions y_n=\sin(n\pi x/L) are already known to be real and mutually orthogonal on [0,L] with weight 1: this theorem is what makes that a special case of one general fact, not a coincidence specific to sines.

The whole argument rests on the orthogonality of \{1, \cos(n\omega_0 t), \sin(n\omega_0 t)\}_{n\ge 1} over one period [-T/2, T/2]: for integers m, n \ge 1, \int_{-T/2}^{T/2}\cos(n\omega_0 t)\cos(m\omega_0 t)\,dt = \begin{cases} T/2 & n=m \\ 0 & n\ne m\end{cases}, \qquad \int_{-T/2}^{T/2}\sin(n\omega_0 t)\sin(m\omega_0 t)\,dt = \begin{cases} T/2 & n=m \\ 0 & n\ne m\end{cases}, and \int_{-T/2}^{T/2}\cos(n\omega_0 t)\sin(m\omega_0 t)\,dt = 0 for every n, m (an odd function integrated over a symmetric interval). Each of these follows from the product-to-sum identities, which turn every integrand into a sum of \cos or \sin of (n\pm m)\omega_0 t, and a full period of any nonzero-frequency sine or cosine integrates to 0.

Fix m \ge 1 and multiply both sides of the series by \cos(m\omega_0 t), then integrate over [-T/2, T/2], exchanging the sum and the integral (justified by the assumed term-by-term integrability): \int_{-T/2}^{T/2} f(t)\cos(m\omega_0 t)\,dt = \frac{a_0}{2}\int_{-T/2}^{T/2}\cos(m\omega_0 t)\,dt + \sum_{n=1}^\infty\Big(a_n\!\int_{-T/2}^{T/2}\!\cos(n\omega_0 t)\cos(m\omega_0 t)\,dt + b_n\!\int_{-T/2}^{T/2}\!\sin(n\omega_0 t)\cos(m\omega_0 t)\,dt\Big). The constant term vanishes (a full period of a nonzero-frequency cosine integrates to 0), every sine-cosine cross term vanishes by orthogonality, and every cosine-cosine term vanishes except n=m, which contributes a_m \cdot T/2. So the right side collapses to a_m T/2, giving a_m = \frac{2}{T}\int_{-T/2}^{T/2} f(t)\cos(m\omega_0 t)\,dt. The same argument with \sin(m\omega_0 t) in place of \cos(m\omega_0 t) isolates b_m the same way, since now every cosine-sine cross term and the constant term vanish, leaving b_m T/2.

For a_0, multiply by 1 instead and integrate: every n \ge 1 term vanishes (a full period of \cos or \sin), leaving \int_{-T/2}^{T/2} f(t)\,dt = \frac{a_0}{2}\cdot T, i.e. a_0 = \frac{2}{T}\int_{-T/2}^{T/2} f(t)\,dt — exactly the general a_n formula evaluated at n=0, since \cos(0) = 1, which is why the constant term is conventionally written a_0/2 rather than a_0.

Only the indicial equation itself is derived here — the three existence cases above are a substantially longer argument (constructing the second solution via reduction of order applied to the first, term by term) that belongs to a dedicated ODE-theory treatment; this course states the three-way split as a fact its own worked examples verify case by case.

Since xp(x) and x^2q(x) are analytic at x=0, write xp(x)=p_0+p_1x+p_2x^2+\cdots and x^2q(x)=q_0+q_1x+q_2x^2+\cdots, so p(x)=p_0/x+p_1+p_2x+\cdots and q(x)=q_0/x^2+q_1/x+q_2+\cdots. Differentiating the trial solution term by term, y=\sum_{n=0}^{\infty}a_nx^{n+r}, \quad y'=\sum_{n=0}^{\infty}(n+r)a_nx^{n+r-1}, \quad y''=\sum_{n=0}^{\infty}(n+r)(n+r-1)a_nx^{n+r-2}. Substitute into y''+p(x)y'+q(x)y=0 and multiply the whole equation by x^2 (clearing every term down to a power of x no worse than x^{n+r}, since x^2\cdot x^{n+r-2}=x^{n+r}, x^2p(x)\cdot x^{n+r-1}=(xp(x))\cdot x^{n+r}, and x^2q(x)\cdot x^{n+r} is already a power series in x^{n+r}). Collecting the coefficient of the lowest surviving power, x^{0+r}=x^r (the n=0 term, using only a_0, p_0=\lim xp(x), and q_0=\lim x^2q(x) since every other contribution to that power carries a strictly positive power of x from the series for p or q): \big[r(r-1)+p_0r+q_0\big]a_0=0. The trial solution requires a_0\neq0 (otherwise r could be raised by one and the same series restarted), so the bracket itself must vanish — exactly the indicial equation above. Every higher power x^{n+r}, n\ge1, produces a recurrence relation expressing a_n in terms of a_0,\dots,a_{n-1} once r is fixed at one of the indicial equation's two roots — the mechanism the three existence cases above are built on.