I recently encountered the expression below for which I was interested in solving for $x$:
\begin{equation} \psi(x+n+1) - \psi(x+1) =y \end{equation}
$\psi$ is the digamma function, $n$ is a positive integer and $x,y>0$. Initially, I solved this numerically but I was interested to see if there was the possibility to solve this, either in closed form or as a series. I tried to find some way by using the relationship to harmonic numbers, e.g. rewriting as
\begin{equation} \psi(x+n+1) - \psi(x+1) = \sum_{k=1}^{n} \frac{1}{k+x} = \int_0^1 \frac{t^x(1-t^n)}{1-t}dt \end{equation}
to see if there was something which I could do, but I didn't notice anything apparent. The only thing I've come up with is that for large values of $x$, $y$ is small and using $\psi(x) \sim\ln(x)$, that $x \approx \dfrac{n}{e^y - 1}$.
I'm happy for any suggestions on how to proceed!
$$\psi(x+n+1) - \psi(x+1) =y$$
For large values of $x$, the lhs write $$\frac{n}{x}-\frac{n (n+1)}{2 x^2}+\frac{n (n+1) (2 n+1)}{6 x^3}-\frac{n^2 (n+1)^2}{4 x^4}+O\left(\frac{1}{x^5}\right)$$ Using series reversion $$x_{(1)}=\frac{n}{y}-\frac{n+1}{2}+\frac{\left(n^2-1\right)}{12 n}y+O\left(y^3\right)$$
Comparing with your approximation $$x_{(2)}= \dfrac{n}{e^y - 1}=\frac{n}{y}-\frac{n}{2}+\frac{n }{12}y+O\left(y^3\right)$$ It is really good since $$x_{(2)}-x_{(1)}=-\frac{1}{2}-\frac{y}{12 n}+O\left(y^3\right)$$
Using the above for $n=1000$ and $y=1$, the exact solution is $581.477$ while $x_{(2)}=582.833$ and $\frac{1000}{e-1}=581.977$.
Edit
We could do better using more terms in the series expansion, then series reversion and transformation of the result into a $[n,n+1]$ Padé approximant.
This would give $$x_{(3)}=\frac{60 (n-1) n^3 -24 n^2 \left(n^2-4\right)y+3 n \left(n^3-3 n^2-9 n-13\right)y^2} {60 (n-1) n^2 y+6 n \left(n^2+11\right)y^2+(n-1) \left(n^2+11\right)y^3 }$$
For the worked example, this gives $x_{(3)}=581.590$.