Solution of $\dfrac{d [\epsilon_i]}{dt} = (\beta \mathbf{A} -\delta \mathbf{I})\left[\epsilon_i\right] - \alpha \mathbf{B} \left[\epsilon_i^2\right]$ I have a problem that can be described by the following equation:
\begin{equation}
\dfrac{d \left[\epsilon_i\right]}{dt} = \left( \beta \mathbf{A} - \delta  \mathbf{I} \right) \left[\epsilon_i\right] - \alpha \mathbf{B} \left[\epsilon_i^2\right],
\end{equation}
where $\left[\mathbf{\epsilon_i}\right]$ is a vector whose components are $\epsilon_i$, $\left[\mathbf{\epsilon_i^2}\right]$ is a vector whose components are $\epsilon_i^2$ and $\mathbf{A}$ and $\mathbf{B}$ are constant matrices with real positive entries. This problem allows the trivial solution $\epsilon_i = 0$, but I am interested on a positive non-trivial solution, i.e. $\epsilon_i > 0$ for some $i$'s (in words, not all elements are zero). In other words, I am interested in finding the parameters $\alpha$, $\beta$ and $\delta$ that allow me to have a positive non-trivial solution. 
For instance, if $\alpha = 0$ I have a standard eigenvalue problem. Note that, considering a real square matrix with positive entries, thus this problem reduces to $\frac{\beta}{\delta} > \lambda_{max}^{-1}(\mathbf{A})$, than such solution exists, where we used the Perron-Frobenius theorem to guarantee that the eigenvalue components are positive. I am interested in a solution where $\alpha > 0$. What would be the constraints involving the relationship between my matrices and my parameters that would allow me to have such solution.
I tried to analyse the system in the steady-state, i.e. $\dfrac{d \left[\epsilon_i\right]}{dt} \rightarrow 0$, but I do not succeed. My first effort was trying to use the Algebraic Riccati equation, but I am not sure whether this would lead to a solution.
Any ideas on how I could tackle this problem?
 A: The case of $2\times 2$ matrices was studied in Matrix Bernoulli Equations. II (2008). It cannot in general be solved by quadratures unless the matrices have a special form. Solvability conditions for higher-dimensional matrices were considered in part I of this series of papers (which I have only found in Russian).
A: If you are interested in conditions guaranteeing that all solutions starting at a point in the non-negative orthant $\mathbb{R}^n_{+}$ remain there, the natural one is the following:  At $(\epsilon_1,\ldots,\epsilon_n)$ the $j$-th coordinate, $F_j$, of the vector field $F$ should be $\ge 0$ for those $j$ for which $\epsilon_j = 0$.  In the linear case (that is, if you take $\alpha = 0$) this translates into ${\beta} A$ being a Metzler matrix, that is, a matrix with all off-diagonal entries non-negative.  Incidentally, the sign and/or the magnitude of $\delta$ are of no relevance here (they add only to the Perron-Frobenius eigenvalue). 
The situation changes dramatically if $\alpha > 0$ (under the assumption that $B$ is a Metzler matrix).  Consider $(\epsilon_1, 0, \ldots, 0)$ with $\epsilon_1 > 0$ (and rather large). We have
$$
F_2(\epsilon_1, 0, \ldots, 0) = {\beta}\, a_{21} \, \epsilon_1 - {\alpha} \, b_{21} \, (\epsilon_1)^2,
$$
which goes to $-\infty$ as $\epsilon_1 \to \infty$.  In particular, for $\epsilon_1$ sufficiently large, 
$$
\left. \frac{d \epsilon_2}{dt}\right\rvert_{t = 0} < 0.
$$
On the other hand, for $(\epsilon_1,\ldots,\epsilon_n)$ close to the origin the linear part is dominant, so the solutions starting in the orthant close to the origin should remain in the orthant, provided that we can show that they remain close to the origin. 
