$$ y'=\pmatrix {3 & 1 \\ -5 & -1}y +e^t \pmatrix {-2 \\ 3}$$ How do I find a solution to the system? I have problem solving the system due to the eigenvalues are complex.This is how far I get
$$x(t)=c_1e^{t+it}+(-2+i)c_1e^{t+it}+c_2e^{t-it}+(-2-i)c_2e^{t-it}$$
You’re almost there for the solution to the homogeneous equation. Recall the superposition principle: every linear combination of solutions to the homogeneous equation is also a solution. Now, you’ve got a pair of independent solutions of the form $e^{\lambda t}\mathbf v$ corresponding to the two complex eigenvalue/eigenvector pairs that you computed. Notice that these two functions are complex conjugates. (You can always arrange for this to be true: If $A$ is a real matrix and $A\mathbf v=\lambda\mathbf v$, then $\overline{A\mathbf v} = A\overline{\mathbf v} = \overline\lambda \overline{\mathbf v}$, so $\overline{\mathbf v}$ is also an eigenvector of $A$ with eigenvalue $\overline\lambda$.) The real and imaginary parts of a complex number are linear combinations of the number and its complex conjugate, so the real and imaginary parts of either of your two primitive solutions also satisfy the homogeneous equation. I’ll leave it to you to verify that these new real-valued solutions to the equation are linearly independent. Thus, a real-valued general solution to the homogeneous equation is $$c_1e^t\begin{bmatrix}\cos t \\ -2\cos t-\sin t\end{bmatrix} + c_2e^t\begin{bmatrix}\sin t \\ \cos t-2\sin t \end{bmatrix}.$$
You could instead take real and imaginary parts earlier in the process. A $2\times2$ real matrix $A$ with complex eigenvalues $\alpha\pm i\beta$ is similar to a matrix of the form $$C = \begin{bmatrix}\alpha&-\beta\\\beta&\alpha\end{bmatrix}$$ with exponential $$e^{tC} = e^{\alpha t}\begin{bmatrix} \cos{\beta t} & -\sin{\beta t} \\ \sin{\beta t} & \cos{\beta t} \end{bmatrix}.$$ We then have $$e^{tA} = e^{\alpha t} P \begin{bmatrix} \cos{\beta t} & -\sin{\beta t} \\ \sin{\beta t} & \cos{\beta t} \end{bmatrix} P^{-1},$$ where the columns of $P$ are the real and imaginary parts of one of the complex eigenvectors. (For a way to compute the exponential of a real $2\times2$ matrix with complex eigenvalues directly from its eigenvalues without computing eigenvectors at all, see this answer.)
As to finding a particular solution of the inhomogeneous equation $\mathbf y'(t) = A\mathbf y(t)+\mathbf b(t)$, a somewhat brute-force approach is to proceed as you might for the one-dimensional equation by making the ansatz $\mathbf y(t) = \exp(tA)\mathbf w(t)$ and seeing where that takes you. Omitting the gory details, this leads to the particular solution $$\mathbf y(t) = \int_0^t \exp[(t-s)A] \,\mathbf b(s) \,ds.$$ This has the handy property that $\mathbf y(0) = 0$, so the general solution to the inhomogeneous equation is obtained by adding the general homogeneous solution to this. For your problem, the components of $\mathbf b(t)$ are themselves exponentials, so evaluating this integral, while tedious, shouldn’t present any serious difficulties.
