Why does variation of parameters work? I was looking at the variation of parameters method, and to be sincere, when I took my differential equations course I felt like too much of it was Hocus Pocus.
For example, there is this
https://tutorial.math.lamar.edu/Classes/DE/VariationofParameters.aspx
I don't know where does (3) comes out from.
It says:
"Now, there is no reason ahead of time to believe that this can be done. However, we will see that this will work out."
My differential equations professor didn't explain it. It wasn't outrageous so I used, simply as a recipe. I would like to have some intuition, some idea behind it.
When you solve some equations in which you need to integrate from both sides, for example
$$
\int f(t) dt = \int g(t) dt
$$
Then you have $F(t)+c_{1}=G(t)+c_{2}$, but you only write $F(t)=G(t)+c$ because $c_{1}$ and $c_{2}$ don't matter, only their difference. I guess something similar is going on there, there are a lots of ways to write that with different functions and you're not interested on the $u_{1}(t)$ and the $u_{2}(t)$ but on a relationship between them, so you choose them accordingly in order to simplify the equation as the relationship will hold.
But I would love to hear some idea, some intuition behind it, as each time that someone ask me for help with differential equations I want to tell them that I know nothing about them, because that's how I really feel.
Thank you for your time.
 A: An interesting way of understanding the variation of parameters is as a multidimensional integrating factor method. You can rewrite your linear system of differential equations as first-order linear differential equation with multiple variables, i.e. as
\begin{equation}
x^\prime(t) + A(t)x(t) = f(t),
\end{equation}
where $x(t),f(t)$ are $n$-dimensional vectors and $A(t)$ is an $n\times n$ matrix (in your case $x(t) = (g(t),g^\prime(t),...,)$ if $g$ is the solution to the higher order linear system). In fact this problem is more general.
Now let $B(t)$ be a fundamental matrix for the corresponding homogeneous system, i.e. the columns of $B(t)$ are $n$ linearly independent solutions to the corresponding equation with $f = 0$. This means that $B(t)$ is invertible for all $t$ as long as the coefficient matrix $A$ is continuous (Abel's theorem), and satisfies
\begin{equation}
B^\prime(t) + A(t)B(t) = 0.
\end{equation}
We can calculate $(B^{-1})^\prime(t)$ by differentiating the identity $B^{-1}(t)B(t) = I$ to get
\begin{equation}
(B^{-1})^\prime(t)B(t) + B^{-1}(t)B^\prime(t) = 0.
\end{equation}
Using that $B^\prime(t) = -A(t)B(t)$, we get
\begin{equation}
(B^{-1})^\prime(t)B(t) = B^{-1}(t)A(t)B(t),
\end{equation}
so that (cancelling the $B(t)$ which is invertible)
\begin{equation}
(B^{-1})^\prime(t) = B^{-1}(t)A(t).
\end{equation}
Thus $B^{-1}(t)$ is an integrating factor for the original equation and we get
\begin{equation}
(B^{-1}x)^\prime(t) = B^{-1}(t)(x^\prime(t) + A(t)x(t)) = B^{-1}(t)f(t).
\end{equation}
Solving for $x(t)$, we get
\begin{equation}
x(t) = B(t)\left(\int_0^t B^{-1}(t)f(t) + C\right),
\end{equation}
where the vector $C$ depends on the initial conditions.
This explains why the solution can be written as a linear combination (with time-varying coefficients) of the fundamental homogeneous solutions (the columns of $B(t)$).
