This is similar but not identical to standard examples in e.g. Paul's Notes, and while the math seems straightforward the results I get disagree with what I get from numerical simulation. Given a 2D system
$$x' = Ax + b$$
my solution is
$$x = c_1e^{\lambda_1t}\eta_1 + c_2e^{\lambda_2t}\eta_2 + A^{-1}b$$
where $\lambda_i$ and $\eta_i$ are eigenvalues and eigenvectors of $A$, respectively. The first two terms in the RHS of the above are straight out of the textbook and indeed my simulation is fine when $b = 0$. Is the $A^{-1}b$ term incorrect?
I realized any system with a constant can be transformed into one without by adding a dimension. $A$ becomes
$$B = \begin{matrix} a_{1,1} & a_{1,2} & b_0 \\ a_{2,1} & a_{2,2} & b_1 \\ 0 & 0 & 0 \end{matrix}$$
and define $x^*$ so that
$$x^*_1(0) = x_1(0) \\ x^*_2(0) = x_2(0) \\ x^*_3(0) = 1$$
then the new system
$$ {x^*}' = Bx^* $$
has the standard solution.