Variation of parameters is a universal method for finding a particular solution to y'' + p(x)y' + q(x)y = f(x). Assume y_p = u₁(x)y₁ + u₂(x)y₂ where y₁, y₂ solve the homogeneous equation, and solve for u₁, u₂ using the Wronskian. Though more computational than undetermined coefficients, this method works for any continuous f(x), making it the universal tool when other methods fail.
To find the general solution of the nonhomogeneous ODE y'' + p(x)y' + q(x)y = f(x), you need a particular solution y_p to add to the general homogeneous solution y_h = c₁y₁ + c₂y₂. The method of undetermined coefficients finds y_p by guessing its form — it works beautifully when f(x) is a polynomial, exponential, sine, or cosine, but fails for f(x) = sec(x), f(x) = ln(x), or any other function that doesn't generate a finite family of derivatives. Variation of parameters solves this by making no assumptions about f(x) at all.
The idea is to promote the constants c₁ and c₂ in the homogeneous solution to functions u₁(x) and u₂(x). You guess y_p = u₁(x)y₁(x) + u₂(x)y₂(x) — the same linear combination that solved the homogeneous equation, but with variable coefficients instead of constants. Substituting into the ODE and imposing a natural simplifying constraint (u₁'y₁ + u₂'y₂ = 0, which eliminates the y'' cross-terms and simplifies the algebra) reduces the problem to a linear system: u₁'y₁ + u₂'y₂ = 0 and u₁'y₁' + u₂'y₂' = f(x). This is a 2×2 system in u₁' and u₂', and your prerequisite Wronskian W = y₁y₂' − y₂y₁' appears directly in the denominator of the solution by Cramer's rule: u₁' = −y₂f/W and u₂' = y₁f/W.
The Wronskian plays a decisive role here. From your prerequisite, you know that W ≠ 0 precisely when y₁ and y₂ are linearly independent — and that is exactly when the 2×2 system has a unique solution. This is why you need a fundamental solution set (two independent solutions) before beginning: without them, the Wronskian would vanish and the method would break down. Once you have u₁' and u₂', you integrate each: u₁(x) = ∫ −y₂f/W dx and u₂(x) = ∫ y₁f/W dx. The integration step is where the actual work lies, and it requires integration by parts when f(x) is complicated.
The method has no restrictions on f(x) beyond continuity on the interval of interest. This is its great advantage over undetermined coefficients, and also explains why it's more computational: undetermined coefficients exploits special structure in f(x) to shortcut the integration; variation of parameters makes no such use of structure and integrates directly. In practice, try undetermined coefficients first when f(x) allows it, and reach for variation of parameters when it does not. The two methods are complementary, not competing.