For any square linear system $\,A\vec x=\vec b\,$ over some field, there exists a unique solution iff $\,\det A\neq 0\,$ , as then we can use the inverse matrix:
$$A\vec x=\vec b\Longleftrightarrow A^{-1}A\vec x=A^{-1}\vec b\Longleftrightarrow A^{-1}\vec b=\vec x $$
As for (a) and your "main question": if $\,\det A=0\,$ one still may have to check whether there are no solutions or infinite solutions (assuming we're working on an infinite field). For example, if the system is homogeneous (over an infinite field) it must have infinite solutions, whereas if the system is non-homogeneous it may have no solutions or several:
$$\begin{cases}x+y=1\\x+y=1\end{cases} \Longleftrightarrow \begin{pmatrix}1&1\\1&1\end{pmatrix}\binom{x}{y}=\binom{1}{1}\longrightarrow\,\,\text{infinite solutions}$$
$$\begin{cases}x+y=1\\x+y=0\end{cases} \Longleftrightarrow \begin{pmatrix}1&1\\1&1\end{pmatrix}\binom{x}{y}=\binom{1}{0}\longrightarrow\,\,\text{no solutions at all}$$
and, of course, in both cases above we have $\,\det A=0\,$
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.
Best Answer
Cramer's rule is useful to solve the linear system $Ax=b$ only if $\text{det}(A)\neq 0$. When $\text{det}(A)\neq 0$ the homogeneous linear system $Ax=0$ has only the zero solution $x=0$. Notice that if the square matrix $A$ is invertible, one has
$$Ax=0 \iff A^{-1}Ax=A^{-1}0 \iff Ix=0 \iff x=0$$ where $I$ is the identity matrix.