I've been wondering for a while, what kind of logic can lead to the method of variation of parameters (i.e. it's trivial to prove, but why would one ever try to find solution this way).
Now I found a convincing interpretation in the book "Обыкновенные Дифференциальные Уравнения" (ODE) by Arnold This interpretation comes from the celestial mechanics and serves a basis for all of the perturbation methods.
In this book the story for the variation of parameters method begins with a sample model: the planets moving around the sun. The first approximation is the motion according to Kepler's laws. Then the objective is to count deviations, caused by planets attracting each other. And the idea is to suppose that planets would still follow Kepler's laws, but the perturbation would make parameters vary over time.
The same idea applies to the inhomogeneous linear DE's. The solution is \(\phi = \phi_h + \phi_p\) the sum of solution for homogeneous system and particular solution for inhomogeneous one.
Let's suppose that \(\phi_h\) is kind of principal (undisturbed) part of the solution, and \(\phi_p\) is some disturbance caused by inhomogenuity. It could happen that disturbance can be described by variations of constants in \(\phi_h\). Now let's just subtitute such solution into equation: voila! That's the intuition that Zadorozhny's course lacked