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
From $x_2'(t)=- x_2(t)\cos t$, you obtain (say, by separation of variables), $x_2(t)=x_2(t_0)e^{\sin t_0}e^{-\sin t}$.
Substitute in the first equation and you obtain $$x_1'(t)+x_1(t)\sin t=x_2(t_0)e^{\sin t_0}te^{-\sin t}$$
The fundamental function of (the homogenous part of) this last ODE is $e^{\cos t}$, and therefore \begin{align}&x_1(t)=x_1(t_0)e^{-\cos t_0}e^{\cos t}+e^{\cos t}\int_{t_0}^t x_2(t_0)e^{\sin t_0}ue^{-\cos u-\sin u}du=\\=&\left(x_1(t_0)e^{-\cos t_0}+x_2(t_0)e^{\sin t_0}\int_{t_0}^0 ue^{-\cos u-\sin u}du\right)e^{\cos t}+x_2(t_0)e^{\sin t_0}e^{\cos t}\int_{0}^t ue^{-\cos u-\sin u}du\end{align}
Which means that $$\begin{pmatrix}x_1(t)\\ x_2(t)\end{pmatrix}=\begin{pmatrix}A e^{\cos t}+Be^{\cos t}\int_0^t ue^{-\cos u-\sin u}\,du\\ Be^{-\sin t}\end{pmatrix}$$
With $A=x_1(t_0)e^{-\cos t_0}+x_2(t_0)e^{\sin t_0}\int_{t_0}^0 ue^{-\cos u-\sin u}\,du$ and $B=x_2(t_0)e^{\sin t_0}$. This means that a possible choice for $\Psi$ is $$\Psi(t)=\begin{pmatrix}e^{\cos t}&e^{\cos t}\int_0^tue^{-\cos u-\sin u}\,du\\ 0&e^{-\sin t}\end{pmatrix}$$
And now you can calculate easily $\Psi(\tau)(\Psi(t))^{-1}$.