Finding the solution to a non-homogeneous matrix exponential.

So you solved the homogeneous system. To now get the solution of the inhomogeneous system you can use the variation of constant technique. Let

$$ \begin{bmatrix}a\\b \end{bmatrix}= e^{Pu}c(u) \longrightarrow \frac{d}{du}\begin{bmatrix}a\\b \end{bmatrix} = Pe^{Pu}c(u) + e^{Pu}c'(u)$$

Thus we have $e^{Pu}c'(u) = \begin{bmatrix}\cos(zu)\\ -\sin(zu)\end{bmatrix}$ and $c(0) = 0$. All that is left to do is to solve this linear system and integrate w.r.t. $u$:

$$c'(u) = e^{-Pu}\begin{bmatrix}\cos(zu)\\ -\sin(zu)\end{bmatrix} = e^{xu}\begin{bmatrix}\cos(yu)& \sin(yu)\\ -\sin(yu) & \cos(yu)\end{bmatrix}^{-1} \cdot\begin{bmatrix}\cos(zu)\\ -\sin(zu)\end{bmatrix} $$

Since the matrix here is orthonormal its inverse is simply its transpose:

$$c'(u) = e^{xu}\begin{bmatrix}\cos(yu)& -\sin(yu)\\ \sin(yu) & \cos(yu)\end{bmatrix} \cdot\begin{bmatrix}\cos(zu)\\ -\sin(zu)\end{bmatrix} = e^{xu}\begin{bmatrix}\cos(yu-zu) \\\sin(yu-zu)\end{bmatrix} $$

Were I used the addition theorems to simplify the expression. Now integrate, using the formulas

\begin{align} \int \sin(\alpha u) e^{\beta u} du &= \frac{\beta\sin(\alpha u) - \alpha\cos(\alpha u)}{\alpha^2+\beta^2}e^{\beta u} +\text{const.}\\ \int \cos(\alpha u) e^{\beta u} du &= \frac{\beta\cos(\alpha u) + \alpha\sin(\alpha u)}{\alpha^2+\beta^2}e^{\beta u} +\text{const.}\\ \end{align}

to get

$$ c(u) = \frac{1}{x^2 + (y-z)^2}e^{xu} \begin{bmatrix} x\cos((y-z)u) + (y-z)\sin((y-z) u)\\ x\sin((y-z) u) - (y-z)\cos((y-z) u) \end{bmatrix} + \text{const.} $$

And to satify the initial condition $c(0)=0$ we need

$$ \text{const.} = -\frac{1}{x^2 + (y-z)^2}\begin{bmatrix} x \\- (y-z) \end{bmatrix}$$

I will use standard notation for differential equations. Your notation is quite unusual and it might lead to unwanted mixing up of formulas.

$$\frac{d}{dt} \begin{bmatrix}x_1\\x_2 \end{bmatrix} = \begin{bmatrix}-a& b\\ -b&-a\end{bmatrix} \begin{bmatrix}x_1\\ x_2\end{bmatrix} + \begin{bmatrix}\cos(\omega t)\\ -\sin(\omega t)\end{bmatrix}$$

In order to determine the solution to the problem, we will compare this to the standard linear control problem. Given a system

$$\dot{\boldsymbol{x}}(t)=\boldsymbol{Ax}(t)+\boldsymbol{Bu}(t)$$

with the initial condition $\boldsymbol{x}(t=t_0)=\boldsymbol{x}_0$.

The solution to this equation is given by

$$\boldsymbol{x}(t)=\exp\left[\left(t-t_0 \right)\boldsymbol{A} \right]\boldsymbol{x}_0+\int_{t=t_0}^{t}\exp\left[(t-\tau)\boldsymbol{A}\right]\boldsymbol{B}\boldsymbol{u}(\tau)d\tau.$$

This formula can be verified by differentiating it and invoking the fundamental theorem of calculus for the integral.

In your case

$$\boldsymbol{A}=\begin{bmatrix}-a& b\\ -b&-a\end{bmatrix},$$ $$\boldsymbol{B}=\begin{bmatrix}1& 0\\ 0 & 1\end{bmatrix},$$ $$\boldsymbol{u}=\begin{bmatrix}\cos\left(\omega t\right)\\-\sin\left(\omega t\right)\end{bmatrix}, \text{ and }$$ $$\boldsymbol{x}_0=\begin{bmatrix}0\\ 0 \end{bmatrix} .$$