Constant term in $\int_{0}^{T}e^{tA}$ equals $-A^{-1}$ I believe the following claims are true, but I cannot prove them. Can anyone provide a proof or a reference?
For any matrix $A$, each entry, $X_{ij}$ of 
$$
X:=\int_{0}^{T}e^{At}dt
$$
can be written as 
$$
\sum_{k=0}^{K_{ij}}a_{ijk}e^{b_{ijk}T}T^{c_{ijk}}
$$
with $K_{ij}<\infty$ and unique, possibly complex,  coefficients $a,b,c$.  If $A$ is invertible, then $A^{-1}=-Y$, where
$$
Y_{ij}:=\sum_{k \text{ s.t. } b_{ikk}=c_{ijk}=0} a_{ijk}.
$$
Note that this is easy to show if $A$ is diagonalizable or if the spectrum of $A$ is strictly contained in the left half plane.
More generally, I believe 
the constant term of 
$$
\int_{0}^{T}e^{At}Ce^{Bt} dt
$$
solves the Sylvester equation
$$
AY+YB=C
$$
whenever a unique solution exists.
 A: The result you give is not exact (see last equation). In fact  if $A$ is a scalar ($1 \times 1$ matrix) :
$$ \text{if} \ A=a, \ \ \int_0^T\exp(ta)dt= T+\tfrac12T^2 a +\cdots \tag{1}$$
And in the general case :
$$\int_0^T\exp(tA)dt= TI_n+\tfrac12T^2 A+\cdots \tag{2}$$
Here is a proof of (2) in the case where $A$ is diagonalizable :
$$A=P\Lambda P^{-1} \ \ \text{where} \ \ \lambda=diag(\lambda_1\cdots \lambda_n)\tag{3}$$
Equivalently :
$$tA=P (t\Lambda) P^{-1} \ \ \text{where} \ \ t\Lambda=diag(t\lambda_1\cdots t\lambda_n)$$
Therefore :
$$\exp(tA)=P M P^{-1} \ \ \text{where} \ \ M=diag(e^{t\lambda_1},\cdots e^{t\lambda_n})$$
$$\int_0^T \exp(tA)dt=P N P^{-1}\tag{4}$$
with 
$$N=diag(\int_0^Te^{t\lambda_1}dt,\cdots \int_0^T e^{t\lambda_n}dt)$$
$$N=diag(\frac{1}{\lambda_1}(e^{\lambda_1T}-1),\cdots \frac{1}{\lambda_n}(e^{\lambda_nT}-1))$$
$$N=\Lambda^{-1}(\exp(T\Lambda)-I_n)\tag{5}$$
Plugging (5) into (4) :
$$\int_0^T\exp(tA)dt=P \Lambda^{-1}(\exp(T\Lambda)-I_n) P^{-1}=P \Lambda^{-1}(T\Lambda+\tfrac12T^2 \Lambda^2+\cdots) P^{-1}$$
giving (2).
A: We change $A,B$ with $-A,-B$.

$\textbf{Proposition}$. Let $A,B,C$ be complex $n\times n$ matrices, $(\lambda_i)_i=spectrum(A),(\mu_i)_i=spectrum(B)$. Assume that , for every $i,j$, $Re(\lambda_i+\mu_j)>0$.

Then the equation $AY+YB=C$ has the unique solution
$X=\int_0^{+\infty}e^{-tA}Ce^{-tB}dt$.
$\textbf{Proof}$. We use the Kronecker product, cf.
https://en.wikipedia.org/wiki/Kronecker_product
Here, we stack the matrices into vectors, row by row.
i) Notice that $f=A\otimes I+I\otimes B^T:Y\mapsto AY+YB$ is one to one. Indeed, $0\notin spectrum(f)=(\lambda_i+\mu_j)_{i,j}$; let $X=f^{-1}(C)$.
Moreover, $-f$ is asymptotically stable (its eigenvalues have a negative real part).
ii) Let $g_t=e^{-tA}\otimes e^{-tB^T}:C\mapsto e^{-tA}Ce^{-tB}$; then $g_t=e^{-tf}$ and
$\int_0^{\infty} g_t(C)dt=(\int_0^{\infty}g_tdt)(C)=f^{-1}(C)=X$. $\square$
