How exactly is the Wronskian tied with Cramer's rule in a system of ODEs? 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.
 A: 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.
