Trick of selecting fundamental matrix in nonhomogeneous differential equations to simplify calculations Problem:
I came across the following differential equation set:
$$
\begin{cases}
 \dot{x}&=3x-2y+t\\
 \dot{y}&=4x-y+t^2\\
\end{cases}
$$
My solution
I followed normal practices and found the general solution to the homogeneous equation $$
\left[ \begin{array}{c}
 x\\
 y\\
\end{array} \right] '=\left[ \begin{matrix}
 3&  -2\\
 4&  -1\\
\end{matrix} \right] \cdot \left[ \begin{array}{c}
 x\\
 y\\
\end{array} \right] \Rightarrow \left[ \begin{array}{c}
 x\\
 y\\
\end{array} \right] =\left[ \begin{array}{c}
 e^t\left( c_1\cos 2t+c_2\sin 2t \right)\\
 e^t\left( \left( c_1-c_2 \right) \cos 2t+\left( c_1+c_2 \right) \sin 2t \right)\\
\end{array} \right] 
$$
I selected $(c_1, c_2)=(1,0)$ and $(0,1)$ to form a fundamental matrix $X=\left[ \begin{matrix}
 e^t\cos 2t&  e^t\sin 2t\\
 e^t\left( \cos 2t+\sin 2t \right)&  e^t\left( \sin 2t-\cos 2t \right)\\
\end{matrix} \right] $, and plugged it into the variation of constants formula to get a particular solution $y(t)=X(t)\int_{0}^{t}\left(X^{-1}(\tau)\cdot \left[ \begin{array}{c}
 \tau\\
 \tau ^2\\
\end{array} \right]\right)\mathrm{d}\tau $, and after a one-page long calculation and the help of Wolframalpha, I got the following solution of $y$:

Question
I wonder what's wrong with my selection since the what wolfram gave me if I typed the original equation set would be like this:

So I guess there must be something in my selection of the fundamental matrix that makes the calculations really complex, but how can I improve it?
Also, what's your favourite way of solving this kind of nonhomogeneous linear differential equation sets?

We want to solve for $x,y$ such that
$$
     \frac{d}{dt}\left[\begin{array}{c}x \\ y\end{array}\right]
    = \left[\begin{array}{cc}3 & -2 \\ 4 & -1\end{array}\right]\left[\begin{array}{c}x \\ y\end{array}\right] +\left[\begin{array}{c}t \\ t^2\end{array}\right].
$$
Let $C$ be the constant coefficient matrix on the right. Then the above may be written as
$$
              \frac{d}{dt}\left[\begin{array}{c}x\\y\end{array}\right]-C\left[\begin{array}{c}x\\y\end{array}\right]=\left[\begin{array}{c} t \\ t^2\end{array}\right] \\
            \frac{d}{dt}\left(e^{-tC}\left[\begin{array}{c} x \\ y\end{array}\right]\right) = e^{-tC}\left[\begin{array}{c} t \\ t^2\end{array}\right]
$$
Integrating over $[0,t]$ and rearranging terms gives
$$
        e^{-tC}\left[\begin{array}{c} x \\ y\end{array}\right] = \left[\begin{array}{c}x_0\\ y_0\end{array}\right]+\int_0^t e^{-sC}\left[\begin{array}{c} s \\ s^2 \end{array}\right]ds \\
        \left[\begin{array}{c} x \\ y\end{array}\right]= e^{tC}\left[\begin{array}{c}x_0 \\ y_0\end{array}\right]+\int_0^t e^{(t-s)C}\left[\begin{array}{c}s \\ s^2\end{array}\right]ds
$$
Now the task is reduced to finding $e^{uC}$ for a real parameter $u$. You can use complex variables to write
$$
     e^{uC} = \frac{1}{2\pi i}\oint_{\gamma} \frac{e^{uz}}{zI-C}dz,
$$
where the contour $\gamma$ is a positively-oriented, simple, closed, rectifiable curve enclosing all eigenvalues of $C$ in its interior. The coefficient matrix $C$ has characteristic polynomial
$$
     p(z)= (z-3)(z+1)+8=z^2-2z+8=(z-4)(z+2)
$$
Therefore, $p(C)=0$, which gives
\begin{align}
    p(z)I&=p(z)I-p(C)\\
               &=(z^2 I-C^2)-2(z I-C) \\
               &=(z I-C)(z I+C)-2(z I-C) \\
               &=(z I-C)(z I+C-2I)
\end{align}
Hence, $z I-C$ is invertible for all $z$ for which $p(z)\ne 0$, and the inverse of $z I-C$ is given by
$$
       (z I-C)^{-1}=\frac{1}{p(z)}((z-2)I+C).
$$
And that gives an integral which may be directly computed by residues at $z=-2,4$, assuming that the contour $C$ encloses $z=2$ and $z=4$ in its interior:
$$
    e^{uC} = \frac{1}{2\pi i}\oint_{\gamma}\frac{e^{uz}}{(z-4)(z+2)}((z-2)I+C)dz.
$$
This is a relatively straightforward integral to evaluate by residues because the integrand has only poles of order $1$ at $z=-2,4$. Hence, $e^{uC}$ is a linear combination of matrices $I$ and $C$:
\begin{align}
   e^{uC}&= \frac{e^{-2u}}{-6}(-4I+C)+\frac{e^{4u}}{6}(2I+C) \\
     &= \left(\frac{2}{3}e^{-2u}+\frac{1}{3}e^{4u}\right)I+\left(-\frac{1}{6}e^{-2u}+\frac{1}{6}e^{4u}\right)C.
\end{align}
Therefore, using the above, the solution is known at this point, even though the fully expanded expression is messy:
$$
        \left[\begin{array}{c} x \\ y\end{array}\right]= e^{tC}\left[\begin{array}{c}x_0 \\ y_0\end{array}\right]+\int_0^t e^{(t-s)C}\left[\begin{array}{c}s \\ s^2\end{array}\right]ds
$$
