Let $A$ be a (finite) Hurwitz matrix.

In this related question of mine, (see also https://en.wikipedia.org/wiki/Lyapunov_equation) it is shown that $$ \int_0^\infty \sum_{j,k} (e^{At})_{ij} Q_{jk} (e^{At})_{\ell,k}dt = W_{i,\ell}, $$ where the matrix $W$ is the solution of some matrix equation $AW+WA^T+Q = 0$. Numerically, the solution of such an equation can be computed in $O(d^3)$, where $d$ is the dimension of my system.

I am wondering if this result can be easily extended to multi-linear maps? In particular, is it possible to compute (numerically) the quantity $$ \int_0^\infty \sum_{j_1,j_2,j_3}Q_{j_1,j_2,j_3}(e^{At})_{i_1,j_1} (e^{At})_{i_2,j_2}(e^{At})_{i_3,j_3}dt $$ by an efficient procedure?

Also, if we go to higher order, like $$ \int_0^\infty \sum_{j_1,j_2,\dots,j_k}Q_{j_1,\dots,j_k}(e^{At})_{i_1,j_1} \dots (e^{At})_{i_k,j_k}dt, $$ can we compute the solution of such an equation with a (relatively efficient) numerical procedure? If yes, how does its complexity grow with $k$?

[EDIT] : after my answer below, there still remains a question whether a more efficient algorithm exists (like the Lyapunov equation for $k=2$).