Liouville's Formula; I think this problem can be solved using Liouville's Formula, like this:
For the sake of brevity in $\LaTeX$, let us agree to denote the coefficient matrix in this problem by $A(t)$, thus:
$A(t) = \begin{bmatrix} \cos^4(t) && 0 && -\sin(2t) \\ \sin(4t) && \sin(t) && -4 \\ -\sin(5t) && 0 && -\cos(t) \end{bmatrix}; \tag{1}$
then if we set
$\mathbf r(t) = \begin{pmatrix} x(t)\\ y(t)\\ z(t) \end{pmatrix}, \tag{2}$
the differential equation in question becomes
$\dot {\mathbf r} (t) = A(t) \mathbf r (t); \tag{3}$
we consider in the usual manner a fundamental solution matrix $X(t, t_0)$ for (3); that is, $X(t, t_0)$ is a $3 \times 3$ matrix function of $t$ satisfying
$\dot X(t, t_0) = A(t) X(t, t_0) \tag{4}$
with
$X(t_0, t_0) = I, \tag{5}$
the $3 \times 3$ identity matrix. It will be noted that the columns of $X(t, t_0)$ are themselves solutions of (3), and that if
$\mathbf r_0 = \begin{pmatrix} x_0 \\ y_0 \\ z_0 \end{pmatrix}, \tag{6}$
then $X(t, t_0)\mathbf r_0$ is the unique solution to (3) with
$\mathbf r(t_0) = \mathbf r_0, \tag{7}$
since
$\dfrac{d}{dt} (X(t, t_0)\mathbf r_0) = \dot X(t, t_0) \mathbf r_0 = (A(t)X(t, t_0)) \mathbf r_0 = A(t)(X(t, t_0) \mathbf r_0); \tag{8}$
it follows that the matrix $X(t, t_0)$ encodes all essential information about the solution space of (3).
We consider the matrix $X(t + 2\pi, t_0)$; we have
$\dot X(t + 2\pi, t_0) = A(t + 2\pi)X(t + 2\pi, t_0) = A(t)X(t + 2\pi, t_0), \tag{9}$
since $A(t)$ is periodic of period $2\pi$: $A(t + 2\pi) = A(t)$; thus $X(t + 2\pi, t_0)$ satisfies the same differential equation as does $X(t, t_0)$, but with initial condition $X(t_0 + 2\pi, t_0)$. Since
$X(t_0 + 2\pi, t_0) = IX(t_0 + 2\pi) = X(t_0, t_0)X(t_0 + 2\pi, t_0), \tag{10}$
it follows from the linearity of (3), and the uniqueness of solutions that for any $t \in \Bbb R$ we have
$X(t + 2\pi, t_0) = X(t, t_0)X(t_0 + 2\pi, t_0) ; \tag{11}$
we see from (11) that further
$X(t + 4\pi, t_0) = X(t, t_0)X(t_0 + 2\pi, t_0)X(t + 2\pi, t_0) = X(t, t_0)X^2(t + 2\pi, t_0) , \tag{12}$
and from here a very simple induction, the completion of which is left to the reader, establishes
$X(t + 2n\pi, t_0) = X(t, t_0)X^n(t_0 + 2\pi, t_0) , \tag{13}$
for any positive $n \in \Bbb Z$. (13) indicates that the long-term growth of the solution matrix $X(t, t_0)$ is intimately tied in with the expansive/contractive properties of the matrix $X(t_0 + 2\pi, t_0)$. In particular, if $X(t_0 + 2\pi, t_0)$ has an eigenvalue $\lambda$ with $\vert \lambda \vert > 1$, and corresponding eigenvector $\mathbf v$, that is
$X(t_0 + 2\pi, t_0)\mathbf v = \lambda \mathbf v \tag{14}$
we find
$X(t + 2n\pi, t_0) \mathbf v = X(t, t_0) X^n(t_0 + 2\pi, t_0)\mathbf v = X(t, t_0)\lambda^n \mathbf v = \lambda^n X(t, t_0) \mathbf v, \tag{15}$
whence, since $\vert \lambda \vert > 1$, the solution $X(t, t_0)\mathbf v$ grows without bound as $t \to \infty$. We further observe that $X(t, t_0)$ is nonsingular for all $t$, and $[t_0, t_0 + 2\pi]$ is compact, so $\Vert X(t, t_0) \mathbf v \Vert$ is bounded below away from $0$ on $[t_0, t_0 + 2\pi]$ by some real $\mu > 0$:
$\Vert X(t, t_0) \mathbf v \Vert > \mu > 0, \; t_0 \le t \le t_0 + 2\pi; \tag{16}$
thus we affirm that $\Vert X(t + 2n\pi, t_0)\mathbf v \Vert > \vert \lambda^n \vert \mu$, which shows that $\Vert X(t + 2n\pi, t_0)\mathbf v \Vert \to \infty$ as $n \to \infty$ independently of $t_0$ or $t$; $X(t, t_0)$ is nonsingular since $X(t_0, t_0) = I$, and the columns of $I$ are linearly independent, and linear independence or dependence of solutions is preserved under the flow of a first-order linear ordinary differential equation.
The preceding discussion indicates that computation of the eigenvalues of $X(t_0 + 2\pi, t_0)$ may be decisive in determining the stability or instability of solutions to (3); however, it is in general a difficult task to explicitly solve (4) for $X(t, t_0)$, and hence equally challenging to find its eigenvalues. But in some cases, such as the present one, progress may be made via a less than direct route; such is it situation here.
If it were possible to evaluate $\det X(t_0 + 2\pi, t_0)$ and to show that $\det X(t_0 + 2\pi, t_0) > 1$, then we could affim that $\vert \lambda \vert > 1$ for at least one eigenvalue of $X(t_0 + 2 \pi, t_0)$ and hence conclude that the system (3) is unstable. Fortunately, for the present problem this is the case, thanks to Liouville's Formula.
Liouville's Formula asserts that, given a system such as (4), $\det X(t_0 + 2\pi, t_0)$ evolves according to the scalar differential equation
$\dfrac{d \det X(t, t_0)}{dt} = \operatorname{Tr}(A) \det X(t, t_0), \tag{17}$
which has an immediate solution, given that $X(t_0, t_0) = I$,
$\det X(t, t_0) = \exp(\displaystyle \int_{t_0}^t \operatorname{Tr}(A(s))ds) \det X(t_0, t_0) = \exp(\displaystyle \int_{t_0}^t \operatorname{Tr}(A(s))ds); \tag{18}$
we have
$\operatorname{Tr}(A(t)) = \cos^4 t + \sin t - \cos t, \tag{19}$
whence
$\displaystyle \int_{t_0}^{t_0 + 2\pi} \operatorname{Tr}(A(s))ds = \int_{t_0}^{t_0 + 2\pi}\cos^4 s ds + \int_{t_0}^{t_0 + 2\pi}\sin s ds - \int_{t_0}^{t_0 + 2\pi}\cos s ds; \tag{20}$
the last two integrals on the right of (20) vanish and we are left with
$\displaystyle \int_{t_0}^{t_0 + 2\pi}\operatorname{Tr}(A(s))ds = \int_{t_0}^{t_0 + 2\pi}\cos^4 s ds; \tag{21}$
we now refer to this MSE post where it is shown that
$\cos^4 t = \dfrac{3 + 4 \cos(2t) + \cos(4t)}{8}, \tag{22}$
whence
$\displaystyle \int_{t_0}^{t_0 + 2\pi} \operatorname{Tr}(A(s))ds = \int_{t_0}^{t_0 + 2\pi} \dfrac{3}{8} ds + \dfrac{1}{2}\int_{t_0}^{t_0 + 2\pi}\cos(2s)ds + \dfrac{1}{2}\int_{t_0}^{t_0 + 2\pi}\cos(4s)ds$
$= \displaystyle \int_{t_0}^{t_0 + 2\pi} \dfrac{3}{8} ds = \dfrac{3\pi}{4}; \tag{23}$
then by (18) we find
$\det X(t_0 + 2\pi, t_0) = e^{3\pi / 4} > 1; \tag{24}$
it now follows from (24) that $X(t_0 + 2\pi, t_0)$ has an eigenvalue $\lambda$ with $\lambda > 1$; hence the system (3) is in fact unstable, i.e., has at least one unbounded solution.
