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.
This is a Euler-Cauchy differential equation. Multiply by $t^2$ to obtain:
$$t^2y’’+4ty’+2y=t^2u.$$
You can determine the solutions of the homogeneous part/fundamental solution by the ansatz $y_h=t^n$.
$$t^2n(n-1)t^{n-2}+4tnt^{n-1}+2t^{n}=0 \quad \text{; if } t\neq 0\implies n^2+3n+2=0 $$
$$\implies (n+2)(n+1)=0 \implies n_1=-1 \qquad n_2 = -2.$$
After using the ansatz and determining the two possible values for $n_{1,2}$ you can assemble the fundamental system
$$x_1=y=c_1t^{-1}+c_2t^{-2} \implies \boldsymbol{X}_{\text{scalar}}=\begin{bmatrix}t^{-1} & t^{-2} \end{bmatrix}$$
$$x_2=\dot{y}=-c_1t^{-2}-2c_2t^{-3}$$
$$\implies \boldsymbol{x}(t)=c_1\begin{bmatrix}t^{-1} \\-t^{-2} \end{bmatrix}+c_2\begin{bmatrix}t^{-2}\\-2t^{-3} \end{bmatrix}$$
The fundamental solution (columns are vectors associated with the constants of integration $c_1$ and $c_2$) is given by:
$$\boldsymbol{X}(t)=\begin{bmatrix}t^{-1}& t^{-2}\\ t^{-2}&-2t^{-3}\end{bmatrix} \implies \boldsymbol{X}^{-1}(t)= \dfrac{1}{\det\boldsymbol{X}(t)}\begin{bmatrix}-2t^{-3}& -t^{-2}\\ -t^{-2}&t^{-1}\end{bmatrix}.$$
An alternative approach for determining the fundamental system of the matrix equation would be to use the scalar fundamental system $\boldsymbol{X}_{\text{scalar}}$ and the following relationship.
$$\boldsymbol{X}(t)=\begin{bmatrix}\boldsymbol{X}_{\text{scalar}}\\\dot{\boldsymbol{X}}_{\text{scalar}}\end{bmatrix} .$$
The state transition matrix $\boldsymbol{\Phi}(t,\tau)$ is given by:
$$\boldsymbol{\Phi}(t,\tau)=\boldsymbol{X}(t)\boldsymbol{X}^{-1}(\tau).$$
An alternative method to directly obtain the state transition matrix is the Peano-Baker series
$${\boldsymbol {\Phi }}(t,\tau )={\boldsymbol{I}}+\int _{\tau }^{t}{\boldsymbol {A}}(\sigma _{1})\,d\sigma _{1}+\int _{\tau }^{t}{\boldsymbol {A}}(\sigma _{1})\int _{\tau }^{{\sigma _{1}}}{\boldsymbol {A}}(\sigma _{2})\,d\sigma _{2}\,d\sigma _{1}$$
$$+\int_{\tau }^{t}{\boldsymbol {A}}(\sigma _{1})\int _{\tau }^{{\sigma _{1}}}{\boldsymbol {A}}(\sigma _{2})\int _{\tau }^{{\sigma _{2}}}{\boldsymbol {A}}(\sigma _{3})\,d\sigma _{3}\,d\sigma _{2}\,d\sigma _{1}+...$$
The general solution is then given by:
$${\boldsymbol {x}}(t)={\boldsymbol {\Phi }}(t,t_{0}){\boldsymbol {x}}(t_{0})+\int _{{t_{0}}}^{t}{\boldsymbol {\Phi }}(t,\tau ){\boldsymbol {B}}(\tau ){\boldsymbol {u}}(\tau )d\tau .$$
As Kwin van der Veen already mentioned your $\boldsymbol{B}(t)$ should be different. In your case it is:
$$\boldsymbol{B}(t)=\begin{bmatrix}0 \\1\end{bmatrix}.$$
Best Answer
Take a simple system $\dot{x} = x$, then $A=1$ and $\phi(t,t_0) = e^{t-t_0}$. That is $A$ is constant, but the system response is time varying. In general this will be the case unless $A=0$.
If you have $\phi(t,t_0)$ in general you have $\dot{\phi}(t,t_0) = A(t) \phi(t,t_0)$, or $A(t) = \phi(t,t_0)^{-1} \dot{\phi}(t,t_0) = \phi(t_0,t) \dot{\phi}(t,t_0)$, so you can determine if $A$ is time varying or not.
Addendum: It is not too hard to show that the system is time invariant iff $\phi(t,t_0) = \phi(t-t_0,0)$ for all $t,t_0$. A little more work shows that in this case $A$ is essentially constant.