Solving Linear Differential Equations: Method of Variation of Constants Imposes Seemingly Arbitrary Constraints

A general way of solving an $$n$$th degree linear (non-homogeneous) differential equation $$L(y)=f$$, given a fundamental set of solutions, is the so-called method of "variation of constants", concisely explained here. Basically, a particular solution $$y_p(x)$$ is proposed as being a linear combination of the fundamental set $$\{y_i(x)\;|\;i\in1 . . n\}$$, but with the constants themselves also a function of $$x$$:

$$y_p(x)=\sum_{i=1}^nu_i(x)y_i(x)$$

... which would then be derived $$n$$ times to be substituted into $$L(y)=f$$.

I first came across this method in my university's textbook on differential equation analysis. At a certain point along the path to the resulting expression for this method (i.e. an integral involving Green's function), the textbook posits the constraints that sums of the following form are assumed to be $$0$$: $$0=\sum_{i=1}^nu_i'(x)y^{(j)}$$ with $$j\in0 . . n-2$$, arguing this "simplifies the derivation work, prevents higher-order derivatives of the constants from entering any equations, and provides the extra equations required for this problem to be solved algebraically".

Now, these arguments all seem backwards to me: in mathematics, one doesn't just "magically" summon equations to reach a desired result, or make the work easier. For example, the equation $$a^2 + b^2 = c^2$$, is easily solved for $$a$$ if I "magically" assume $$b=4$$ and $$c=5$$; three unknowns, three linearly independent equations.

I utilised various sources to check if these constraints were general practice, and indeed, they are, but there doesn't seem to be a concise explanation for them:

• My professor at university (she wrote the textbook) told me she doesn't know why we do it; "it just works".
• The Wikipedia page linked earlier argues "For this purpose, one adds the constraints ...", but doesn't explain why this would be allowed.
• The Wikipedia page on the variation of constants argues "(...) where the $$u_i(x)$$ are differentiable functions which are assumed to satisfy the conditions ...", but doesn't explain why this would be allowed.
• The authors of the book this Stack Exchange user was reading on the matter, "suggest to set ... for the sake of convenience".
• The respondent in that same Stack Exchange thread argues, for the case of $$n=2$$, that "If we let the terms of the first derivative of $$u_1$$, and $$u_2$$ to subsist, then $$y^{(2)}$$ will have terms of the second derivative of $$u_1$$ and $$u_2$$, and we will achieve nothing in terms of simplicity. What was a second order ODE with one unknown ($$y$$), will become a second order ODE with two unknowns ($$u_1$$ and $$u_2$$), with the lack of a second equation. This is why, in the construction, one asks for ... If this condition is fulfilled, we achieve one very important thing: to reduce the non-homogeneous second order ODE, to a system of two first order ODEs, which, incidentally, we know how to solve exactly (modulus integration). This would be one rational explanation in choosing the above condition."
• The respondent in this Stack Exchange thread argues, in a comment, "It's what we call an "ansatz", a guess based on experience, knowledge, that we know will help our proof."

This all seems to be very non-rigorous mathematics: "we need extra constraints, so we create a equations that incidentally make all kinds of convenient results pop out of the original problem".

Surely, there must be a more rigorous justification as to why these constraints may be used, and why these constraints are used?