This answer is an attempt to provide a solution for a hypothetical user that does not know the general formula beforehand and is unaware of any packages. The answer below considers that the user is aware of the following: - If $f$ is an expandable function that can be considered on the real line (for example $f(x)=(1+x)^{-1}$ or $f(x)=\exp(x)$), then, for a matrix variable $A$ and a real number valued variable $a$, $f(A)$=`Replace`[$f(a),a \rightarrow A$]. This point allows us to use Series for real valued variables. - For $A$ and $B$ two matrices $\left(A.B\right)^{-1}=B^{-1}.A^{-1}$. This point is used to place the initial expression into the form $f(x)$ where $f(x)=(1+x)^{-1}$ and $ x=A^{-1}.B $ Then if we define an identity tensor for the dot product: Dot[s___, id, g___] ^:= Dot[s, g] The formula may be obtained with the following code that attempts to do simple replacements sequentially (var below is a formal wrapper that formally replaces a matrix with a real valued variable) : Fold[ReplaceAll, Inverse[A + B] , { Inverse[w_] :> Inverse[A . TensorExpand[Inverse@A . w]] , Inverse[a_ . b_] :> Inverse[b] . Inverse[a], MatrixPower[_, 0] -> id, Inverse[id + s_] :> ReplaceAll[ Normal@Series[(1 + epsilon*var[s])^(-1), {epsilon, 0, 3}] /. epsilon -> 1 , { var[a_]^m_:>Dot@@ConstantArray[a,m] , var->Identity , 1 -> id } ] } ]//TensorExpand out: `(* -MatrixPower[A, -1] . B . MatrixPower[A, -1] + MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] - MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] + MatrixPower[A, -1] *)`