Suppose we have the following system:
\begin{equation} \dot{x} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -1 & 0 & -2 & 0 \\ \end{bmatrix}x + Bu \end{equation}
\begin{equation} y = Cx \end{equation} I want to find conditions on $B \in \mathbb{R}^{4\text{x}1}$ in order this system's zero state response to be bounded for all $C \in \mathbb{C}^{1\text{x}4}$ for $u(t) = constant$. But isn't this the definition of BIBO stability? If it is, then since the eigenvalues of this system are $\lambda = \{-j,j\}$ this system cannot be BIBO stable since there will be no poles with negative real parts. If this is not equivalent to BIBO stability, then how can $y(t) = Ce^{At}Bu$ be bounded? There are $2$ jordan blocks with sizes $2$ and it would be almost impossible to get rid of $te^{jt}$ terms.
First all, if you want the output to be bounded for all $C\in\mathbb{C}^{1\times 4}$, then this is equivalent to saying that it is bounded for all $C$ in $\{e_1^T,\ldots,e_4^T\}$ where the $e_i$'s are the vector of the natural basis for $\mathbb{C}^4$.
We have that
$$\exp(At)=\left(\begin{array}{cccc} \cos\left(t\right)+\frac{t\,\sin\left(t\right)}{2} & \frac{3\,\sin\left(t\right)}{2}-\frac{t\,\cos\left(t\right)}{2} & \frac{t\,\sin\left(t\right)}{2} & \frac{\sin\left(t\right)}{2}-\frac{t\,\cos\left(t\right)}{2}\\ \frac{t\,\cos\left(t\right)}{2}-\frac{\sin\left(t\right)}{2} & \cos\left(t\right)+\frac{t\,\sin\left(t\right)}{2} & \frac{\sin\left(t\right)}{2}+\frac{t\,\cos\left(t\right)}{2} & \frac{t\,\sin\left(t\right)}{2}\\ -\frac{t\,\sin\left(t\right)}{2} & \frac{t\,\cos\left(t\right)}{2}-\frac{\sin\left(t\right)}{2} & \cos\left(t\right)-\frac{t\,\sin\left(t\right)}{2} & \frac{\sin\left(t\right)}{2}+\frac{t\,\cos\left(t\right)}{2}\\ -\frac{\sin\left(t\right)}{2}-\frac{t\,\cos\left(t\right)}{2} & -\frac{t\,\sin\left(t\right)}{2} & -\frac{3\,\sin\left(t\right)}{2}-\frac{t\,\cos\left(t\right)}{2} & \cos\left(t\right)-\frac{t\,\sin\left(t\right)}{2} \end{array}\right).$$
One can observe that the only way this remains bounded is by choosing $B$ such that the terms in $t$ cancel each other on each row. This is the case if and only
$$B=\begin{bmatrix}\alpha\\\beta\\-\alpha\\-\beta\end{bmatrix}$$
for any $\alpha,\beta\in\mathbb{R}$. We have that
$$\exp(At)B = \begin{bmatrix} \alpha\cos(t) + \beta\sin(t)\\ \beta\cos(t) - \alpha\sin(t)\\ -\alpha\cos(t) - \beta\sin(t)\\ \alpha\sin(t) - \beta\cos(t) \end{bmatrix}$$
and we clearly have that $C\exp(At)B$ is bounded for all $C\in\mathbb{C}^{1\times 4}$. However, boundedness is not enough, what is important here is that the integral of this expression is bounded for all $t\ge0$ since we will be looking at the step response. This is obviously the case here.
Note that this is not equivalent to asking whether the system is BIBO stable because we are restricting ourselves to the case of constants inputs. The output will not be bounded for the input $u(t)=\sin(t)$ for instance even if $B$ is chosen as above. The concept of stability considered here is, in fact, weaker than BIBO stability.
Additional questions from the comments:
Is it possible to find the same result without having to compute the exponential $\exp(At)$?
Yes, this is definitely possible since $A$ and $\exp(At)$ have the same eigenvectors. The idea is that the matrix $B$ should be in the span of the eigenvectors of the eigenvalues on the imaginary axis but should not be in the span of the generalized eigenvectors.
In other words let $u_i$ be those eigenvectors, and it is necessary and sufficient that $B=\mathrm{span}\{u_1,\ldots\}$.
In the present case, we have that $(A-jI)u_1=0$ implies that $u_1=(j,-1,-j,1)$ and that $(A+jI)u_2=0$ implies that $u_2=(-j,-1,j,1)$. If we want to restrict to real matrices, we just choose $B$ to be in the span of $(u_1+u_2)/2$ and $(u_1-u_2)/(2j)$, and we find exactly the same result as above.
This easily generalizes to a general matrix $A$ where $B$ is in the span of