Conditions for bounded solutions to inhomogeneous system of linear ODEs Consider the simple autonomous, non-homogeneous system of linear ODE's
$$\dot{\mathbf{x}}(t) = M\mathbf{x}(t) + \mathbf{b}$$
where $\mathbf{x} :\mathbb{R}\to\mathbb{C}^n$ is a vector valued function of time, $M$ is a constant $n\times n$ complex matrix, and $\mathbf{b}\in\mathbb{C}^n$ is constant.
My main question is the following:

Q: What are the necessary conditions placed on $M$ and $\mathbf{b}$ such that all solutions $\mathbf{x}(t)$ are bounded, in norm, for all $t$? That is, all solutions satisfy $\|\mathbf{x}(t)\|\leq R$ for all $t$, with  $R\in(0,\infty)$ fixed.

When $\mathbf{b} = 0$, I know the answer: the eigenvalues $\lambda$ of $M$ must satisfy $\text{Re}(\lambda)\leq 0$ and the Jordan blocks (in a Jordan normal form of $M$) corresponding to any eigenvalue with $\text{Re}(\lambda)=0$ must size $1\times 1$.
For $\mathbf{b} \neq 0$, can the Jordan normal form of $M$ alone determine whether all solutions are bounded? Or, does $\mathbf{b}$ play a non-trivial role in diagnosing this notion of "stability"? I think I can prove that all solutions are bounded whenever $M$ is diagonalizable and has eigenvalues that lie strictly in the left half plane ($\text{Re}(\lambda)<0$) but am not sure how to move beyond.
I should also say that I am aware of other notions of stability such as Lyapunov stability, asymptotic stability, etc, and I know various theorems relating the Jordan normal form of $M$ to these notions. I'm simply having trouble finding results relating to this boundedness criterion I state.
Thanks to anyone who may be able to answer or point me in the right direction!
 A: The solution can be calculated as
$$ \dot {\mathbf{x}}(t) - M{\mathbf x}(t) = \mathbf b   $$
$$ (e^{-tM}{\mathbf x}(t))' = e^{-tM}\mathbf b $$
$$ \mathbf x(t) = e^{tM}\mathbf x_0 + \int_0^t e^{(t-s)M}\mathbf b\ ds $$
It is well known that $M = S D S^{-1}$, where $D$ is block diagonal consisting of blocks
$$ J=\begin{pmatrix} \lambda & 1 & 0 &\ldots &0 \\ 0 & \lambda & 1 & \ldots \\  0 & 0 & \lambda & \ldots    \\
 &  & & \ddots & \\ 0&0 &0& \ldots & \lambda \end{pmatrix}, $$
furthermore
$$ \exp(tJ) = \begin{pmatrix} e^{t\lambda} & t\lambda e^{t\lambda} & \ldots &\frac{(t\lambda)^{m-1}}{(m-1)!} e^{t\lambda} \\ 0 & e^{t\lambda}& \ldots &\frac{(t\lambda)^{m-2}}{(m-2)!}e^{t\lambda} \\
& & \ddots & \\
0 & 0 & \ldots & e^{t\lambda} \end{pmatrix} .$$
Now, letting $\mathbf{y}(t) = S^{-1} \mathbf{x}(t) $ we have
$$ \mathbf{y}'(t) = S^{-1}\mathbf{x}'(t) = S^{-1} M \mathbf{x}(t) + S^{-1}\mathbf b = D\mathbf{y}(t) + \widetilde{\mathbf b}. $$
For an arbitrary index $k$ we have
$$ y_k(t) = e^{t\lambda_k}\sum_{j=0}^{l_k-1} \frac{(t\lambda_k)^{j}}{(j)!} y_{0k} +\int_0^t e^{(t-s)\lambda_k}\sum_{j=0}^{l_k-1} \frac{((t-s)\lambda_k)^{j}}{(j)!} \widetilde b_k\ ds  $$
if $\lambda_k \neq 0$. In case $\lambda_k = 0$
$$ y_k(t) = \sum_{j=0}^{l_k-1} \frac{(t\lambda_k)^{j}}{(j)!} y_{0k} +\int_0^t \sum_{j=0}^{l_k-1} \frac{((t-s)\lambda_k)^{j}}{(j)!} \widetilde b_k\ ds . $$
For fixed $j$ we have
$$ \int_0^te^{\lambda_k(t-s)}(\lambda_k(t-s))^j\ ds = t^{j+1}\lambda_k^j\int_0^1 e^{\lambda_k t u}u^j\ du = \big(\lambda_k^{j-1} t^j +O(t^{j-1})\big)e^{\lambda_kt} $$
meaning that the stability of the solution of fully determined by the eigenstructure of $M$.
