I want to solve the system $$\begin{cases} x' = 4x - 13y + 2e^{2t}\cos(3t) \\ y' = x \end{cases}$$ using undetermined coefficients. For the associated homogeneous solution, I get $$ \vec{x_h} = Ae^{2t}\begin{bmatrix}13\cos(3t) \\ 2\cos(3t) + 3\sin(3t)\end{bmatrix} + Be^{2t}\begin{bmatrix}13\sin(3t) \\ 2\sin(3t) - 3\cos(3t)\end{bmatrix}. $$
My first guess for the particular solution was $$ \vec{x_p} = \begin{bmatrix}a_1 \\ a_2\end{bmatrix}e^{2t}\cos(3t) + \begin{bmatrix}b_1 \\ b_2\end{bmatrix}e^{2t}\sin(3t). $$ This was a risky guess, since it closely resembles parts of $\vec{x_h}$, but I don't know the exact conditions under which the default guess fails when working with systems. Sure enough, I ended up with an inconsistent system for $a_1$, $a_2$, $b_1$, and $b_2$.
I looked at the math I used to solve the corresponding second-order equation, $y'' - 4y' + 13y = 2e^{2t}\cos(3t)$, and saw that the guess $y_p = ate^{2t}\cos(3t) + bte^{2t}\sin(3t)$ was successful, so my next guess was $$ \vec{x_p} = \begin{bmatrix}a_1 \\ a_2\end{bmatrix}te^{2t}\cos(3t) + \begin{bmatrix}b_1 \\ b_2\end{bmatrix}te^{2t}\sin(3t). $$ I used the exponential shift rule (which I'm assuming/hoping is valid for systems?) to get $$ e^{2t}(D + 2)\left[\begin{bmatrix}a_1 \\ a_2\end{bmatrix}t\cos(3t) + \begin{bmatrix}b_1 \\ b_2\end{bmatrix}t\sin(3t)\right] = \begin{bmatrix}4 & -13 \\ 1 & 0\end{bmatrix}\left(\begin{bmatrix}a_1 \\ a_2\end{bmatrix}te^{2t}\cos(3t) + \begin{bmatrix}b_1 \\ b_2\end{bmatrix}te^{2t}\sin(3t)\right) + \begin{bmatrix}2e^{2t}\cos(3t) \\ 0\end{bmatrix},$$ but this once again led to an inconsistent algebraic system.
Edit: The system I get here is $$\begin{cases}3b_1 + 2a_1 = 4a_1 - 13a_2 \\ -3a_1 + 2b_1 = 4b_1 - 13b_2 \\ a_1 = 2 \\ b_1 = 0 \\ 3b_2 +2a_2 = a_1 \\ -3a_2 + 2b_2 = b_1 \\ a_2 = 0 \\ b_2 = 0\end{cases}$$
After computing the solution with Wolfram, it appears that the required guess was the sum of the two above, $$ \vec{x_p} = \begin{bmatrix}a_1 \\ a_2\end{bmatrix}te^{2t}\cos(3t) + \begin{bmatrix}b_1 \\ b_2\end{bmatrix}te^{2t}\sin(3t) + \begin{bmatrix}c_1 \\ c_2\end{bmatrix}e^{2t}\cos(3t) + \begin{bmatrix}d_1 \\ d_2\end{bmatrix}e^{2t}\sin(3t), $$ but this seems to break the general pattern. For example, if the forcing term in a second-order linear ODE is $e^t$, and $y_h = e^t$ is an associated homogeneous solution, then the guess $y_p = ate^t$ will work; you don't need $y_p = ate^t + be^t$. The corresponding second-order equation serves as another example. According to this answer, the guess for systems where all but one component of the forcing term are $0$ should be the same as if that non-zero component were the forcing term of a single linear ODE.
Why does this problem, by contrast, require the original terms to be kept when the new terms with an extra factor of $t$ are introduced? Additionally, why does the corresponding second-order equation, $y'' - 4y' + 13y = 2e^{2t}\cos(3t)$, provide misleading information about what the guess should be?