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


In [this related question of mine](https://mathoverflow.net/questions/255474/integral-of-the-entrywise-square-of-the-exponential-of-a-matrix), (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$).



Higher order Lyapunov equation
Inspired by Will's answer of [Integral of the entrywise square of the exponential of a matrix](https://mathoverflow.net/questions/255474/integral-of-the-entrywise-square-of-the-exponential-of-a-matrix) I realized that the above equation can also be computed by using the Kronecker product.


Let $$M= A \oplus A \oplus A = A \otimes I\_n \otimes I\_n + I\_n \otimes A \otimes I\_n + I\_n \otimes I\_n \otimes A.$$
In terms of indices, this means: 
$$M\_{i\_1i\_2i\_3,j\_1j\_2j\_3} = A\_{i\_1j\_1}\delta\_{i\_2j\_2}\delta\_{i\_3j\_3} + 
\delta\_{i\_1j\_1}A\_{i\_2j\_2}\delta\_{i\_3j\_3}+ \delta\_{i\_1j\_1}\delta\_{i\_2j\_2}A\_{i\_3j\_3}$$
We have: 
$$\exp(Mt) = \exp(At)\otimes \exp(At) \otimes \exp(At).$$
In terms of indices this means that: 
$$\exp(Mt)\_{i\_1i\_2i\_3,j\_1j\_2j\_3} = \exp(At)\_{i\_1,j\_1}\exp(At)\_{i\_2,j\_2}\exp(At)\_{i\_3,j\_3}$$ 
In particular, 
$$ \int\_0^\infty \exp(Mt)dt = -M^{-1},$$
which gives the solution to the above integral in $O(n^9)$. 


This can be easily generalized to higher order, with a complexity in $O(n^{3k})$. For $k=2$, it is less efficient than solving the Lyapunov equation (which can be done in $O(n^3)$).


