I assume $X\ge0$ means $u^\top X u\ge0$, and that $B$ *is definite positive* $$\inf_{\|u\|=1} u^\top B\, u:=\beta>0.$$ I also assume matrix norms are the Euclidean operator norms.

Compute the differential by the chain rule, as suggested in comments by F.Poloni:
$$Df(X)H=(AX^{-1}A^{\top}+B)^{-1}AX^{-1}\cdot H\cdot X^{-1}A^{\top}(AX^{-1}A^{\top}+B)^{-1}$$
$$=(A^{\top}+XA^{-1}B)^{-1}\cdot H \cdot(A+BA^{-\top}X)^{-1}=$$
$$=\big(B^\top+Y\big)^{-1}B^\top A^{-\top}\cdot H\cdot A^{-1}B^\top\big(B^\top+Z)^{-1}, $$

where $Y:=B^{\top}A^{-\top} X A^{-1}B\ge0 $ and $Z:=BA^{-\top}XA^{-1}B^\top\ge0$, conjugated to $X\ge0$. Thus for any unit norm vector $u\in\mathbb{R}^n $
$$\big\|\big(B^\top+Y\big)u\big\|\ge u^\top\big(B^\top+Y\big)u\ge u^\top B u \ge\beta$$ and
$$\big\|\big(B^\top+Z)u\big\|\ge u^\top\big(B^\top+Z)u \ge u^\top B u \ge\beta.$$
Hence
$$\big\|\big(B^\top+Y\big)^{-1}\big\|\le \beta^{-1}$$ and
$$\big\|\big(B^\top+Z)^{-1}\big\|\le \beta^{-1}.$$
Therefore
$$\|Df\|_\infty\le\|B\|^2\|A^{-1}\|^2\beta^{-2}$$
which is also a Lipschitz constant for $f$, since its domain is convex, $\{X\ge0\}$.

$$*$$

*Rmk.* The above bounds on the $L_2$ operator norms (or others matrix norms) could be improved, but not up to $$\big\|\big(A^{\top}+XA^{-1}B\big)^{-1}\big\|\le\big\|A^{-1}\big\|,$$
even for *symmetric* definite positive matrices.
Take e.g. $n=2$ and
$$A=:I=\begin{bmatrix}
1 & 0 \\
0 & 1
\end{bmatrix}\quad X:=\begin{bmatrix}
5/2 & 1 \\
1 & 1/2
\end{bmatrix}\quad B:=\begin{bmatrix}
1 & -1/2 \\
-1/2 & 1/2
\end{bmatrix}$$

Then
$$\big(A^{\top}+XA^{-1}B\big)^{-1}=\big(I+XB\big)^{-1}= \begin{bmatrix}
4/15 & 4/15\\
-4/15 & 16/15
\end{bmatrix}$$
whose maximum singular value is $\displaystyle{2\over 15}\sqrt{29}+{2\over 5}>1=\|A^{-1}\|$. The same holds for the Frobenius, and other common entry-wise norms (due to the coefficient $16/15>1$).

$v^\top A v >0$ for nonzero vectors $v$? Or$A$ has positive real eigenvalues? $\endgroup$ – Pietro Majer Jun 18 at 11:52it's not easy to find the derivative of a function of a matrix over a matrixActually it's quite easy. All you need here is the formula $(X+H)^{-1} = X^{-1} - X^{-1}HX^{-1} + O(\|H\|^2)$, which follows from the Neumann series. Applying it twice you should get your derivative. $\endgroup$ – Federico Poloni Jun 18 at 13:15