is there are specific way to solve coupled first-order differential equations with coefficients varying? suppose I have "n" coupled differential equation represented by the matrix, Y• =  A Y  , where Y• is the column matrix containing first derivatives, namely, y1•(t), y2•(t), ...  yn•(t) . A is a square matrix whose each element contains some function dependent on "t" (not constants) and Y is the column matrix containing the solution set, namely, y1(t), y2(t), ...  yn(t) .

If A, contained constants, then its easy to solve by Matrix Exponential method or Eigen-Value method. But, if it contains some varying functions, then is there any approach to solve this. Please, direct me to a good reference, if possible, with an example.
 A: If your matrices commute at different times, that is $A(t)A(s) = A(s)A(t)$, then the solution is $Y(t) = \exp \Big( \int_0^t A(s) \mathrm{d}s \Big) Y(0)$. Here $\exp$ is just an ordinary matrix exponential.
If they do not commute, then you might be able to use Dyson formula (Dyson series). In the realm of quantum mechanics you have the following proposition:
Let $\mathscr{H}$ be a Hilbert space, $A: \mathbb{R} \to \mathscr{B}(\mathscr{H})$ a strongly continuous Hermitian operator valued function (i.e. each $A(t)$ is bounded and Hermitian). Then there is a unique solution of
$$ i \partial_t \psi(t) = A(t) \psi(t), \quad \psi(0) \in \mathscr{H}. $$
If $A$ is continuous w.r.t. the operator norm, then the solution is of the form
$$ \psi(t) = T \exp \Big( - i \int_0^t A(\tau) \mathrm{d}\tau \Big) \psi(0), $$
where $T \exp$ is the so called "time-ordered exponential".
I suspect that the assumption of Hermiticity is not important (you need it for unitarity of the solution). Boundedness, on the other hand, is essential. But this is no problem in yout case of a matrix.
