Approach 1: You already have
$$\tag 1 x(t) = c_1e^{-t} \begin{bmatrix}1\\0\end{bmatrix} + c_2e^{-2t} \begin{bmatrix}1\\1\end{bmatrix}, x[0] = \begin{bmatrix}3\\1\end{bmatrix}$$
Substitute in $t = 0$, equate terms to the IC and find $c_1 = 2, c_2=1$.
Approach 2: Fundamental State Transition matrix
From $(1)$, we can write:
$$\Phi(t) = \begin{bmatrix}e^{-t}& e^{-2t}\\0&e^{-2t}\end{bmatrix}$$
Find:
- $\Phi^{-1}(0)$
- $e^{A t} = \Phi(t) \cdot \Phi^{-1}(0)$
$$x(t) = e^{A t} \cdot x[0]$$
Of course the answers are the same. Also, there are many approaches to solving these.
If $A(t_1)$ and $A(t_2)$ commute for all $t_1$ and $t_2$, so $A(t_1)\,A(t_2) = A(t_2)\,A(t_1)$, then you could use
$$
\Phi(t,0) = e^{\int_0^t A(\tau)\,d\tau} = e^{B(t)}
$$
which indeed for the time-invariant case comes down to $e^{A\,t}$. In order to see why, one can use the Taylor expansion as a way of calculating the matrix exponential
$$
e^{B(t)} = I + B(t) + \frac{1}{2!} B(t)^2 + \frac{1}{3!} B(t)^3 + \cdots
$$
The time derivative of this expression should be equal to $A(t)\,\Phi(t,0)$. Taking the time derivative of this Taylor expansion gives
$$
\frac{d}{dt}e^{B(t)} = \dot{B}(t) + \frac{1}{2!} \left(\dot{B}(t)\,B(t) + B(t)\,\dot{B}(t)\right) + \frac{1}{3!} \left(\dot{B}(t)\,B(t)^2 + B(t)\,\dot{B}(t)\,B(t) + B(t)^2\,\dot{B}(t)\right) + \cdots
$$
when using that $\dot{B}(t) = A(t)$ and that $B(t)$ can be seen as a linear sum of all $A(\tau)\,\forall \tau\in[0,t]$, then if $A(t_1)$ and $A(t_2)$ commute for all $t_1$ and $t_2$ it follows that $A(t)$ and $B(t)$ should commute as well. This allows for $\dot{B}(t) = A(t)$ to be factored out for all term of the derivative of the Taylor series
$$
\frac{d}{dt}e^{B(t)} = A(t)\left(I + \frac{1}{2!} 2\,B(t) + \frac{1}{3!} 3\,B(t)^2 + \cdots\right) = A(t)\left(I + B(t) + \frac{1}{2!} B(t)^2 + \cdots\right)
$$
which is just equal to $A(t)\,e^{B(t)}$.
However your $A(t)$ does not satisfy this condition, so this solution for the state transition matrix can not be used. I am not aware of another way one could for finding an analytical solution for the state transition matrix besides the Peano-Baker series. In this document it is mentioned at page 40 that the state transition matrix is almost always obtained via numerical integration. That document is from 1991. So if back then numerical integration would be preferred over analytical solutions, then nowadays with much more computation power at our disposal it should probably still be the preferred method.
However it can be noted that once you have obtained a solution for $\Phi(t,0)\,\forall\,t\in[0,T]$ with $T$ the period of $A(t)$, then all further solutions can be derived from this using the following expression recursively
$$
\Phi(t+T,0) = \Phi(t,0)\,\Phi(T,0).
$$
From this you can also say something about stability, namely
$$
\Phi(n\,T,0) = \Phi(T,0)^n\quad \forall\ n\in\mathbb{Z}
$$
so if $\Phi(T,0)$ is a Schur matrix (all eigenvalues inside the unit circle) then the system will be stable.
Best Answer
Alternative method 1: A different method is to use the Laplace transform. If you apply the Laplace transform to
$$\dot{\boldsymbol{x}}=\boldsymbol{Ax}$$
you will obtain
$$s\boldsymbol{X}(s)-\boldsymbol{x}(0)=\boldsymbol{AX}(s) \implies \boldsymbol{X}(s)=\left[s\boldsymbol{I}-\boldsymbol{A} \right]^{-1}\boldsymbol{x}(0).$$
If we apply the Laplace inverse we obtain
$$\boldsymbol{x}(t)=\mathcal{L}^{-1}\left\{\left[s\boldsymbol{I}-\boldsymbol{A} \right]^{-1}\right\}\boldsymbol{x}(0).$$
The state transition matrix $\boldsymbol{\Phi}$ is then given by
$$\boldsymbol{\Phi}(t,0)=\mathcal{L}^{-1}\left\{\left[s\boldsymbol{I}-\boldsymbol{A} \right]^{-1}\right\}$$
Alternative method 2: You can solve the system explicitly determining the eigenvalues $\lambda_i$ (Note, that we have three distinct eigenvalues which are on the diagonal of the system matrix) and eigenvectors $\boldsymbol{v}_i$. The general solution is then given as
$$\boldsymbol{x}(t)=c_1\boldsymbol{v}_1\exp(\lambda_1 t)+c_2\boldsymbol{v}_2\exp(\lambda_2 t)+c_3\boldsymbol{v}_3\exp(\lambda_3 t).$$
The fundamental solution $\boldsymbol{Y}$ is given as the matrix
$$\boldsymbol{Y}(t)=\begin{bmatrix} \boldsymbol{v}_1\exp(\lambda_1 t) & \boldsymbol{v}_2\exp(\lambda_2 t) & \boldsymbol{v}_3\exp(\lambda_3 t)\\ \end{bmatrix}.$$
The state Transition matrix $\boldsymbol{\Phi}$ is given as
$$\boldsymbol{\Phi}(t,t'=0)=\boldsymbol{Y}(t)\boldsymbol{Y}^{-1}(t'=0).$$
Note, that we have to invert a $3\times 3$ matrix in the second term.