I have an ODE for an unknown $x(t):[0,\infty)\to\mathbb R^n$ of the following form: $$ x_i'(t)=a_i^\top x(t) + x(t)^\top Q_i x(t), $$ for $i\in\{1,\ldots,n\}$. Here, the vectors $a_i\in\mathbb R^n$ form the columns of an $n\times n$ negative semidefinite matrix $A\in\mathbb R^{n\times n}$, and $Q_i\succeq0$ for all $i\in\{1,\ldots,n\}$. In words, the right-hand side of the ODE is a convex quadratic function of the state. We can write $x'=Ax+\mathcal Q(x)$, where each entry of $\mathcal Q$ is a convex homogeneous quadratic in $x(t)$.
Suppose $x'=Ax$ is well-approximated using backward (implicit) Euler integration to advance $t$, or even that we're willing to do exponential integration for this linear term. Is there a numerical integration method that might be effective in this case for timestepping the nonlinear ODE above?
As a simple example, integration based on variation of parameters could "bake in" solution of the linear part of the ODE. Or perhaps there is a variational version of implicit integration that would lead to a convex problem in our setting.