When when want to find particular solutions to a LDE of order 2 like :
$$ y'' + ay' + by = f(x) $$
I see in a lot of source that says that after finding two solutions for the homogeneous equations $y_1, y_2$, we want to find two functions $\lambda(x), \mu(x)$ that satisfies the system :
$$ \Biggl \{ \begin{gather*}\lambda'(x) y_1 + \mu'(x) y_2 = 0 \\ \lambda'(x) y_1' + \mu'(x) y_2' = f(x) \end{gather*} \qquad (1) $$
And after we can apply Cramer's Rule to find their derivated form, and then integrate them.
But I can't find where does the $(1)$ system comes from.
May it be with the Method of the Variation of Constants ?
The idea behind variation of parameters is to replace the general solution $c_1 y_1 + c_2 y_2$ to the homogeneous system with $c_1(x) y_1 + c_2(x) y_2$. When you plug that into an inhomogeneous equation $y''+a(x) y' + b(x) y = f$, you get
$$c_1'' y_1 + 2 c_1' y_1' + c_1 y_1'' + c_2'' y_2 + 2 c_2' y_2' + c_2 y_2'' + a c_1' y_1 + a c_1 y_1' + a c_2' y_2 + a c_2 y_2' + b c_1 y_1 + b c_2 y_2=f.$$
That's a big mess, but a lot of it can be canceled out by using the fact that both of the $y_i$ satisfy the homogeneous equation. That means $c_i y_i'' + ac_i y_i' + b c_i y_i = 0$ so you have
$$c_1'' y_1 + 2 c_1' y_1' + c_2'' y_2 + 2 c_2' y_2' + a c_1' y_1 + a c_2' y_2 = f.$$
The main reason why this was a good idea in the first place is that the $c_i$ now only appear through their derivatives, so we can do reduction of order. Explicitly, let $d_i=c_i'$ and then you have
$$d_1' y_1 + 2 d_1 y_1' + d_2' y_2 + 2 d_2 y_2' + a d_1 y_1 + a d_2 y_2 = f.$$
Note that given the $d_i$, we can pick the $c_i$ to be any antiderivative of them we want, because shifting $c_i$ by a constant just amounts to moving homogeneous solutions between the particular solution and the homogeneous part of the general solution.
That's one linear first order ODE in two unknown functions, so we can more or less arbitrarily impose another one as long as it is linearly independent of this one. The traditional choice is $d_1 y_1 + d_2 y_2 = 0$. That quantity being a constant means its derivative is zero, so as a bonus we get $d_1' y_1 + d_1 y_1' + d_2' y_2 + d_2 y_2'=0$. Plugging both of those in, the original equation becomes just $d_1 y_1'+d_2 y_2'=f$. This is is linearly independent of our "extra" equation as long as $y_1$ and $y_2$ were linearly independent.
Now since the $y_i$ are assumed known, these are just algebraic equations for $d_i$, so we just solve the algebraic system and then get the $c_i$ by integration.