Symbolic inverse of a linear combination of two matrices This question is very much related to Inverse of the sum of matrices .
Let $d\ge 1$ be an integer and let $a\in(-1,1)$ and $b\in(-1,1)$ be real numbers.
Let ${\bf C}:=\left( f(|i-j|) \right)_{i,j=1}^{d+1}$ be a symmetric positive definite matrix for dimension $(d+1)$ and finally let ${\bf U}:= \left( b^{l-1} a^{k-l-1} 1_{l\ge 1} 1_{k \ge l+1}\right)_{l,k=0}^d$. It appears that the matrix ${\bf U}$ has a rank equal to $(d-1)$.
The problem is now to find:
\begin{equation}
\tilde{{\bf C}}(k):=\left( {\bf C}^{-1} + k \cdot {\bf U} \right)^{-1} = ?
\end{equation}
Ideally we would like to find the function above symbolically as a function of the parameter $k$.
Now, by using the invertible matrix theorem we have shown that:
\begin{eqnarray}
\lim_{k\rightarrow \infty} \tilde{{\bf C}}(k) = 
\left(
\begin{array}{ccccc}
f(0) +\kappa_{1,1} & 0 & \cdots & 0 & f(d) + \kappa_{1,d+1} \\
f(1) +\kappa_{2,1} & 0 & \cdots & 0 & f(d-1) + \kappa_{2,d+1}\\
0 & 0 & \cdots & 0 & 0\\
\vdots \\
0 & 0 & \cdots & 0 & 0
\end{array}
\right)_{(d+1)\times (d+1)}
\end{eqnarray}
where
\begin{eqnarray}
\kappa_{1,1} &=& \frac{1}{{\mathfrak D}} \cdot \sum\limits_{j=0}^{d-2}(-1)^{j+1} f(j+2) \cdot \mbox{det} \left( f(|k-l-1+1_{k\ge 1}+1_{k\ge j+1}|) \right)_{k,l=0}^{d-2}\\
\kappa_{1,d+1} &=& \frac{1}{{\mathfrak D}} \cdot \sum\limits_{j=0}^{d-2}(-1)^{j+1} f(d-2-j) \cdot \mbox{det} \left( f(|k-l-1+1_{k\ge 1}+1_{k\ge j+1}|) \right)_{k,l=0}^{d-2}\\
\kappa_{2,1} &=& \frac{1}{{\mathfrak D}} \cdot \sum\limits_{j=0}^{d-2}(-1)^{j+1} f(j+2) \cdot \mbox{det} \left( f(|k-l+1_{k\ge j+1}|) \right)_{k,l=0}^{d-2}\\
\kappa_{2,d+1} &=& \frac{1}{{\mathfrak D}} \cdot \sum\limits_{j=0}^{d-2}(-1)^{j+1} f(d-2-j) \cdot \mbox{det} \left( f(|k-l+1_{k\ge j+1}|) \right)_{k,l=0}^{d-2}
\end{eqnarray}
and
\begin{equation}
{\mathfrak D}:= \mbox{det} \left( f(|i-j+1|)\right)_{i,j=1}^{d-2}
\end{equation}
Note 1: Our result is consistent with the Woodbury matrix theorem https://en.wikipedia.org/wiki/Woodbury_matrix_identity. Indeed our inverse is a sum of the original matrix ${\bf C}$ and a rank-$(d-1)$ correction to it (at least in the limit $k\rightarrow \infty$) as the theorem states.
Now, it is plausible that finding our function ${\tilde{\bf C}}(k)$ in full symbolic form is a very hard task and in fact is not really necessary because we can always use the theorem in question to compute our function numerically. Yet the temptation is strong to find at least the $1/k$-corrections to the large-$k$ behavior stated above.
How can we do that?
 A: In the following we replace $d+1$ by $d$ without loss of generality.
Since the matrix $C$ is invertible we can write:
\begin{equation}
\left( C^{-1}+ k \cdot U \right)^{-1}= \left(C^{-1} \cdot (1+k \cdot C\cdot  U) \right)^{-1} = \left(1+k \cdot C \cdot U\right)^{-1} \cdot C
\end{equation}
The Cayley-Hamilton theorem for matrix $A:=C \cdot U$ reads:
\begin{equation}
\sum\limits_{l=0}^d {\mathcal a}_{d-l} (-1)^l A^l = 0
\end{equation}
Using this theorem it is not hard to see that :
\begin{equation}
\left(1+k \cdot C \cdot U\right)^{-1} = \sum\limits_{l=0}^{d-1} {\bf O}(k)_{d,d-l} \cdot [C\cdot U]^l
\end{equation}
where 
\begin{eqnarray}
[{\bf O}(k)]^{-1} = {\mathbb 1}_{d\times d} + k \cdot \left( \begin{array}{rrrr} {\mathcal a}_1 & -{\mathcal a}_2 & {\mathcal a}_3 & \cdots & (-1)^{d-2} {\mathcal a}_{d-1} &(-1)^{d-1} {\mathcal a}_d \\ 1& 0 & 0 & \cdots & 0 & 0\\0 & 1 & 0 & \cdots & 0 & 0\\ 0 & 0 & 1 & \cdots & 0 & 0\\
\vdots & \vdots & \vdots & \cdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots &1 & 0\end{array}\right)_{d\times d}
\end{eqnarray}
and the coefficients $\left( {\mathcal a}_l \right)_{l=1}^d$ are the rotational invariants of the matrix $A$, i.e:
\begin{eqnarray}
{\mathcal a}_1 &:=& Tr(A) \\
{\mathcal a}_2 &:=& \frac{1}{2} \left( Tr(A)^2 - Tr(A^2) \right) \\
{\mathcal a}_3 &:=& \frac{1}{6} \left( Tr(A)^3 - 3Tr(A^2)Tr(A)+2 Tr(A^3) \right) \\
{\mathcal a}_4 &:=& \frac{1}{4!} \left( Tr(A)^4 - 6 Tr(A)^2 Tr(A^2) + 8 Tr(A) Tr(A^3) + 3 Tr(A^2)^2-6 Tr(A^4)\right) \\
\vdots\\
{\mathcal a}_{d} &=& \det(A)
\end{eqnarray}
It is not hard to see that the matrix ${\bf O}(k)$ can be found in closed form. One way to do it is to fix $d$, generate that matrix in Mathematica, then symbolically invert and see the pattern. Then the final result is the following:
\begin{eqnarray}
\left(1+k \cdot C \cdot U\right)^{-1} = \sum\limits_{l=0}^{d-1} (-1)^l 
\left(
\frac{\sum\limits_{j=0}^{d-l-1} {\mathcal a}_j k^{l+j}}
{ \sum\limits_{j=0}^d {\mathcal a}_j k^j} \right) \cdot [C\cdot U]^l
\end{eqnarray}
subject to ${\mathcal a}_0=1$.
So far this result is quite general. Now in our case it appears that the matrix $A$ is singular,i.e. both ${\mathcal a}_d = {\mathcal a}_{d-1}=0$.
We take this into account and we calculate the large-$k$ limit of our matrix inverse:
\begin{eqnarray}
\left(C^{-1}+ k \cdot U \right)^{-1} &=& 
k \sum\limits_{l=0}^{d-1} (-1)^l \frac{a_{d-1-l} k^{d-2} + a_{d-l-2}k^{d-3} 1_{l\le d-2} + O(k^{d-4})}{a_{d-2} k^{d-2} + a_{d-3} k^{d-3} + O(k^{d-4})} \cdot [C.U]^l \cdot C \\
&=& 
k \sum\limits_{l=0}^{d-1} (-1)^l \left( \frac{a_{d-l-1}}{a_{d-2}} + \frac{a_{d-l-2} a_{d-2} - a_{d-l-1} a_{d-3}}{a_{d-2}} \cdot \frac{1}{a_{d-3}+a_{d-2} k} + O(\frac{1}{k^2})\right) \cdot [C.U]^l \cdot C \\
& \underset{k\rightarrow \infty}{\rightarrow} & \frac{1}{a_{d-2}} \sum\limits_{l=0}^{d-2} (-1)^l a_{d-l-2} \cdot [C.U]^l\cdot C
\end{eqnarray}
where in the last line we used the fact that $\sum\limits_{l=0}^{d-1} (-1)^l a_{d-l-1} \cdot A^l = 0 $.
Now we have checked numerically that the limit above is exactly equal to the limit stated in the formulation of the question.
