How to obtain matrix from summation inverse equation

This does not seem a simple matrix equation to solve. Computationally, my first attempt would be with Newton's method, even if it takes $O(k^6)$ per iteration, where $k$ is the size of the matrices. The Jacobian of the map in the LHS is $$ L_B f[H] = sum_{i=1}^n (A_i+B)^{-1}H(A_i+B)^{-1}, $$ and to solve the equation $L_Bf[H] = Y$ for $H$ given $Y$ you need to convert it to a $k^2 times k^2$ linear system (there are faster algorithms to solve this linear matrix equation for $n=2$, but I do not think there is anything better otherwise).

If you need to solve it for dimensions for which this is unfeasible, then I would try turning it into a fixed-point equation, for instance $$ B = left(K - sum_{i=2}^n (A_i+B)^{-1}right)^{-1} - A_1, $$ and then hope that the iteration $$ B_{m+1} = left(K - sum_{i=2}^n (A_i+B_m)^{-1}right)^{-1} - A_1, $$ converges.

The scalar version of this equation is a secular equation, but searching for this term I found nothing interesting for matrix arguments.