Let $A$ and $B$ be two positive definite matrices. I would like to find the following integral
$$\int_0^\infty (A x +I)^{-1} B (A x +I)^{-1} dx$$
where $x$ is a scalar and integration is done cell-wise.
Solution when $A$ and $B$ commute
\begin{align} \int_0^\infty (A x +I)^{-1} B (A x +I)^{-1} dx &=\int_0^\infty (A B^{-1}x +B^{-1})^{-1} (A x +I)^{-1} dx \\ &=\int_0^\infty ( B^{-1} A x +B^{-1})^{-1} (A x +I)^{-1} dx\\ &=B \int_0^\infty ( A x +I )^{-1} (A x +I)^{-1} dx\\ &=B A^{-1} \end{align}
The last expression can be inferred from this question.
Now, how to do this for the general case?
At least we can track down its effect under a nicely chosen basis.
I will assume that $A$ and $B$ are $d\times d$ matrices over $\Bbb{C}$ and that $A$ is positive-definite and self-adjoint. (The case with the base field equal to $\Bbb{R}$ is essentially the same.)
Spectral calculus solution. By the spectral theorem, there exists an orthonormal basis $v_1, \cdots, v_d$ in $\Bbb{C}^d$ which consists of eigenvectors of $A$ corresponding to eigenvalues $\lambda_1, \cdots, \lambda_d$, respectively. If $T(B)$ denotes the integral, then with the standard inner product $\langle \cdot, \cdot \rangle$ on $\Bbb{C}^d$, we have
\begin{align*} \langle v_i, T(B) v_j \rangle &= \int_{0}^{\infty} \langle (I + xA)^{-1}v_i, B (I + xA)^{-1}v_j \rangle \, dx \\ &= \int_{0}^{\infty} \langle (I + x\lambda_i)^{-1}v_i, B (I + x\lambda_j)^{-1}v_j \rangle \, dx \\ &= \langle v_i, Bv_j \rangle \int_{0}^{\infty} \frac{dx}{(1+x\lambda_i)(1+x\lambda_j)} \\ &= \mu_{ij}\langle v_i, Bv_j \rangle, \end{align*}
where $\mu_{i,j}$ is given by
$$ \mu_{ij} = \begin{cases} \frac{\log \lambda_i - \log \lambda_j}{\lambda_i - \lambda_j}, & \lambda_i \neq \lambda_j \\ \frac{1}{\lambda_i}, & \lambda_i = \lambda_j \end{cases} $$
A slightly different way of writing this is as follows: let $P_i$ be the orthogonal projection onto the vector $v_i$, i.e., $P_i w = \langle v_i, w \rangle v_i$, or $P_i = |v_i \rangle\langle v_i |$ using the bra-ket notation. Then
$$ T(B) = \sum_{i,j} \mu_{ij} P_i B P_j. $$
Commutator solution. Let us write $R = (I + xA)^{-1}$ and $\mathcal{L}_A B = [A, B] = AB - BA$. Then it is easy to check that
$$ BR = RB + x R(\mathcal{L}_A B)R. $$
Indeed, this follows by writing $x \mathcal{L}_A B = (I+xA)B - B(I+xA)$. Using this relation repeatedly, we get
$$ BR = \sum_{k=0}^{n-1} x^k R^{k+1} \mathcal{L}_{A}^{k} B + x^n R^n (\mathcal{L}_A^{k}B) R, $$
where $\mathcal{L}_A^k$ is the $k$-fold composition of the operator $\mathcal{L}_A$. If $x$ is sufficiently small, we can let $n\to\infty$ to obtain
$$ BR = \sum_{k=0}^{\infty} x^k R^{k+1} \mathcal{L}_{A}^{k} B. $$
Then for $\epsilon$ sufficiently small, we can integrate term by term to obtain
\begin{align*} \int_{0}^{\epsilon} (1+xA)^{-1}B(1+xA)^{-1} \, dx &= \sum_{k=0}^{\infty} \left( \int_{0}^{\epsilon} x^k (1+xA)^{-k-2} \, dx \right) \mathcal{L}_{A}^{k} B \\ &= \sum_{k=0}^{\infty} \frac{1}{k+1} (\epsilon^{-1} + A)^{-k-1} \mathcal{L}_{A}^{k} B. \end{align*}
Now assume that one of the following conditions hold:
Then the series converges uniformly for all $\epsilon > 0$ and thus we can let $\epsilon \to \infty$ to obtain
$$\int_{0}^{\infty} (1+xA)^{-1}B(1+xA)^{-1} \, dx = \sum_{k=0}^{\infty} \frac{1}{k+1} A^{-k-1} \mathcal{L}_{A}^{k} B. $$