Multivariable PVI Let $\beta\in\mathbb{R}$. Solve
$$\left\{\begin{array}{l} y_{1}^{\prime}=y_{1}+\beta y_{2}+12 e^{x},\\ y_{2}^{\prime}=y_{2}+6 e^{x},\\ y_{1}(0)=1,\\ y_{2}(0)=-1. \end{array}\right.$$
Attempt:
The problem can be written as
$$\vec{y}'=\begin{pmatrix}
        1 &  \beta\\
        0 &  1\\
        \end{pmatrix}\vec{y}+\begin{pmatrix}
        12e^{x}\\
        6e^{x}
        \end{pmatrix}.$$
The solution to $\vec{y}'=\begin{pmatrix}1 &  \beta\\0 &  1\\ \end{pmatrix}\vec{y}$ is given by
$$\vec{y}_{c}(x)=e^{xA}\vec{c}.$$
Since
$$A=\begin{pmatrix}
        1 &  \beta\\
        0 &  1\\
        \end{pmatrix}=I_{2}+\beta\begin{pmatrix}
        0 &  1\\
        0 &  0\\
        \end{pmatrix},$$
It follows that
$$\begin{align*} &e^{xA}=e^{xI_{2}}e^{x\beta\begin{pmatrix} 0 &  1\\ 0 &  0\\ \end{pmatrix}}=\begin{pmatrix} e^{x} &  e^{\beta x}\\ 1 &  e^{x}\\ \end{pmatrix}\\ \implies&\vec{y}_{c}(x)=\begin{pmatrix} e^{x} &  e^{\beta x}\\ 1 &  e^{x}\\ \end{pmatrix}\vec{c}. \end{align*}$$
The particular solution is given by
$$\begin{align*}
        \vec{y}_{p}&=e^{xA}\int e^{-xA}\begin{pmatrix}12e^{x}\\6e^{x}\end{pmatrix}dx\\
        &=e^{xA}\int \begin{pmatrix} e^{-x} &  e^{-\beta x}\\ 1 &  e^{-x}\\ \end{pmatrix}\begin{pmatrix}12e^{x}\\6e^{x}\end{pmatrix}dx\\
        &=\begin{pmatrix}
        6(x+2e^{x})e^{\beta x}+6e^{x}\left ( e^{\frac{x-\beta x}{1-\beta}}+2x \right ) \\
        6\left ( e^{\frac{x-\beta x}{1-\beta}}+2x \right )+6e^{x}\left ( x+2e^{x} \right )
        \end{pmatrix}.
    \end{align*}$$
Therefore, the general solution is
$$\vec{y}(x)=e^{xA}\vec{c}+\begin{pmatrix}
        6(x+2e^{x})e^{\beta x}+6e^{x}\left ( e^{\frac{x-\beta x}{1-\beta}}+2x \right ) \\
        6\left ( e^{\frac{x-\beta x}{1-\beta}}+2x \right )+6e^{x}\left ( x+2e^{x} \right )
        \end{pmatrix}.$$
What is wrong with this? When I solve for $\vec{c}$ I get $\begin{pmatrix}-17\\-19\end{pmatrix}$, but when I compute $\vec{y}(0)$, I don't get $\begin{pmatrix}1\\-1\end{pmatrix}$.
 A: You have a defective matrix with a repeated eigenvalue and need a generalized eigenvector. I will skip this approach since you are trying to use the nilpotent matrix.
Using your approach and from @amd's answer at Matrix exponential of non diagonalizable matrix?.
Repeated Eigenvalue: Let $G=A-\lambda I$, where $\lambda$ is the eigenvalue. By the Cayley-Hamilton theorem, $(A-\lambda I)^2=0$, so $G$ is nilpotent. Since $G$ and $\lambda I$ commute, $\exp(tA)=\exp(\lambda tI)\exp(tG)$, but $\exp(tG)=I+tG$ (expand using the power series), so $$\exp(tA) = \mathrm e^{\lambda t}(I+tG).$$
This approach is nice as it needs no eigenvectors. For your problem, we have
$$e^{xA} = e^x\left(I + x \left(
\begin{array}{cc}
 0 & \beta \\
 0 & 0 \\
\end{array}
\right)\right) = e^x \left(
\begin{array}{cc}
 1 & \beta x \\
 0 & 1 \\
\end{array}
\right) $$
You can also do this with eigenvalues and eigenvectors and other ways too. For example, this is a decoupled system and we could have solved the second equation and then substituted that into the first.
You should try the rest of the problem again after this slight correction.
The final answer I get is
$$\begin{align} y_1(x) &= e^x \left(3 \beta x^2-\beta x+12 x+1\right) \\ y_2(x) &= e^x (6x - 1) \end{align}$$
