If we have a nonhomogeneous system of linear differential equations: $$\frac{d\vec{y}}{dx}=A\vec{y}+\vec{b}$$ where: $$\frac{d\vec{y}}{dx}=\begin{bmatrix}y_1'\\y_2'\\y_3'\end{bmatrix},\vec{y}=\begin{bmatrix}y_1\\y_2\\y_3\end{bmatrix},A=\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{bmatrix}, \vec{b}=\begin{bmatrix}b_1(x)\\b_2(x)\\b_3(x)\end{bmatrix}$$ the particular solution involving these three functions could be written: $$y_{particular}=c_1(x)e^{\lambda_1x}\vec{u_1}+c_2(x)e^{\lambda_2x}\vec{u_2}+c_3(x)e^{\lambda_3x}\vec{u_3}$$ where $\lambda_i, \vec{u_i}$ are the eigenvalues and eigenvectors of the constant matrix $A$ (I chose it constant to make things simpler.)
When trying to evaluate these $c_i(x)$ functions, I was told to evaluate them by: $$c_i(x)=\int\frac{W_i}{W}dx$$ Where $W$ is the Wronskian determinant whose columns are the elements of the vectors $c_i(x)e^{\lambda_3i}\vec{u_i}$, and $W_i$ is the Wronskian determinant whose $i$th column was replaced by the contents of the $\vec{b}$ vector.
My questions are: How does the integral show up? Why does the process look like Cramer's rule? How did we lead up to that result? Where did this 'Wronskian' determinant come from?
If I have missed something, please point it out and I'll add it to the post.
Thanks in advance.
Write the general solution of the homogeneous linear ODE $$ \tag{1} \frac{\mathrm{d}\vec{y}}{\mathrm{d}x} = A(x) \vec{y} $$ in the matrix form as $$ \Phi(x; x_0) \vec{c}, $$ where $\vec{c} = \mathrm{col}(c_1, \dots, c_n)$ is a constant matrix (that is, independent of $x$), and $\Phi$ is the state-transition matrix such that $\Phi(x_0, x_0) = I$ (the identity matrix).
Now we are looking for a general solution of the nonhomogeneous linear ODE $$ \tag{2} \frac{\mathrm{d}\vec{y}}{\mathrm{d}x} = A(x) \vec{y} + \vec{b}(x) $$ in the form $$ \Phi(x, x_0) \vec{c}(x). $$ Incidentally, here is the source of the term "variation of constants". We look for conditions for the above function of $x$ to be a solution of $(2)$, that is, to have $$ \frac{\mathrm{d}}{\mathrm{d}x}(\Phi(x, x_0) \vec{c}(x)) = A(x) (\Phi(x, x_0) \vec{c}(x)) + \vec{b}(x). $$ But $$ \frac{\mathrm{d}}{\mathrm{d}x}(\Phi(x, x_0) \vec{c}(x)) = \frac{\mathrm{d}}{\mathrm{d}x} \Phi(x, x_0) \cdot \vec{c}(x) + \Phi(x, x_0) \cdot \vec{c}'(x) $$ and $$ \frac{\mathrm{d}}{\mathrm{d}x} \Phi(x, x_0) = A(x) \Phi(x, x_0), $$ so after cancellation we obtain $$ \Phi(x, x_0) \cdot \vec{c}'(x) = \vec{b}(x), $$ which is, for each fixed $x$, a Cramer (because $\det{\Phi(x, x_0)} = W(x, x_0) \ne 0$) system of $n$ linear equations with $n$ unknowns, $c'_1(x), \dots, c'_n(x)$. The unknowns are given by the Cramer rule. We have thus obtained the derivative of the matrix function $\vec{c}(x)$, so we have to integrate that derivative. And that is the source of the integral.
The "Wrońskian" determinant was named after a nineteenth-century Polish philosopher and mathematician, Józef Maria Hoene-Wroński.
EDIT: I changed slightly the terminology and provided a reference to a Wikipedia entry.