Inverting $\frac{\xi}{2}(1+\tanh(\xi))=\lambda$ using the Lagrange-Burmann Theorem For my quantum mechanics homework, I developed the transcendental equation $\frac{\xi}{2}(1+\tanh(\xi))$ for the well-posedness of symmetric potential formed from two delta functions. The professor encourages us to use a numerical tool to solve the equation
$$\frac{\xi}{2}(1+\tanh(\xi))=\lambda$$
for $\xi(\lambda)$; however, I was curious if the Lagrange inversion theorem could be employed instead.
Taking $B_n$ to represent a Bernoulli number, we have
$$
\begin{eqnarray}
\tanh x &=& x - \frac {x^3} {3} + \frac {2x^5} {15} - \frac {17x^7} {315} + \cdots = \sum_{n=1}^\infty \frac{2^{2n}(2^{2n}-1)B_{2n} x^{2n-1}}{(2n)!}, \left |x \right | < \frac {\pi} {2} \\
\end{eqnarray}
$$
could I pose the transcendental equation as
$$
\frac{\xi}{2} + \sum_{n=1}^\infty \frac{2^{2n-1}(2^{2n}-1)B_{2n} \xi^{2n}}{(2n)!} = \lambda
$$
Excerpting from Wikipedia, if $g$ shall be the inverse of $f$, where $f$ is given by a formal power series as
$$f(w) = \sum_{k=0}^{\infty}f_k\frac{w^k}{k!}$$
$$g(z) = \sum_{k=0}^{\infty}g_k\frac{z^k}{k!}$$
then
$$g_1=\frac{1}{f_1}$$
and
$$g_n = \frac{1}{f_1^n}\sum_{k=1}^{n-1}(-1)^k n^{(k)} \mathcal B_{n-1,k}\left(\frac{f_2}{2f_1},\frac{f_3}{3 f_1},\dots,\frac{f_{n-k+1}}{(n-k+1)f_1}\right)$$
where $n^{(k)}$ is the rising factorial and $\mathcal B$ is a Bell polynomial.
Question: Are there any further simplifications that I can use? Currently, the presence of Bell polynomials seems discouraging as the performance of this algorithm. Also, can I work around the stipulation $|\xi|<\frac{\pi}{2}$? There are some solutions that exist outside of that regime for sufficiently high $\lambda$. Could I partition the function $\frac{\xi}{2}(1+\tanh(\xi))$ into intervals of length $\frac{\pi}{2}$ and apply the result to each of them?
 A: We can re-write the equation as
$$
\lambda  = \frac{\xi }{{1 + {\rm e}^{ - 2\xi } }},
$$
i.e.,
$$
- 2\lambda  = \frac{{ - 2\xi }}{{1 + {\rm e}^{ - 2\xi } }}.
$$
Consequently,
$$
2\lambda {\rm e}^{ - 2\lambda }  = \frac{{2\xi }}{{1 + {\rm e}^{ - 2\xi } }}\exp \left( {\frac{{ - 2\xi }}{{1 + {\rm e}^{ - 2\xi } }}} \right)
$$
or
$$
2\lambda {\rm e}^{ - 2\lambda }  = \frac{{2\xi }}{{{\rm e}^{2\xi }  + 1}}\exp \left( {\frac{{2\xi }}{{{\rm e}^{2\xi }  + 1}}} \right).
$$
Hence, in terms of the Lambert $W$-function,
$$
\frac{{2\xi }}{{{\rm e}^{2\xi }  + 1}}=W(2\lambda {\rm e}^{ - 2\lambda } ).
$$
But from the original equation,
$$
\frac{{2\xi }}{{{\rm e}^{2\xi }  + 1}} = 2(\xi  - \lambda ).
$$
Thus,
$$\boxed{
\xi  = \lambda  + \frac{1}{2}W(2\lambda {\rm e}^{ - 2\lambda } ).}
$$
For example, if $|2\lambda {\rm e}^{ - 2\lambda }|<\frac{1}{\mathrm{e}}$ (which includes $\lambda \ge 0$), then
$$
\xi  = \lambda  + \frac{1}{2}\sum\limits_{n = 1}^\infty  {\frac{{( - n)^{n - 1} }}{{n!}}(2\lambda {\rm e}^{ - 2\lambda } )^n } .
$$
Indeed, this follows by taking the principal brach of $W$ and using its standard Maclaurin series. Alternatively, by $\mathrm{A}038049$,
$$
\xi  =  - \frac{1}{2}\sum\limits_{n = 1}^\infty  {\left( {\sum\limits_{k = 0}^n {\binom{n}{k}k^{n - 1} } } \right)\frac{{( - 2\lambda )^n }}{{n!}}} 
$$
provided $|\lambda|  < \frac{1}{2}W\!\left( {\frac{1}{{\rm e}}} \right) = 0.13923227 \ldots$ ($W$ being the principal branch).
A: Here is a later answer using Lagrange reversion.You can simplify the equation into:
$$y(\tanh(y)+1)=z\implies y=\sum_{n=1}^\infty \frac{z^n}{n!}\left.\frac{d^{n-1}}{dt^{n-1}}(\tanh(t)+1)^{-n}\right|_{t=0},z=2\lambda$$
The result requires finding $n$th derivatives. Rearranging and applying the binomial theorem:
$$\left.\frac{d^{n-1}}{dt^{n-1}}(\tanh(t)+1)^{-n}\right|_{t=0}=2^{-n}\left.\frac{d^{n-1}}{dt^{n-1}}(e^{-2t}+1)^n\right|_{t=0}=2^{-n}\sum_{k=0}^n\binom nk\left.\frac{d^{n-1}}{dt^{n-1}}e^{-2kt}\right|_{t=0}$$
After removing the $n=1$ term for no convergence problems. A simplification and radius of convergence is due to @Gary:
$$\boxed{y(\tanh(y)+1)=z\implies y=z-\frac12 \sum_{n=2}^\infty\sum_{k=0}^n\frac{k^{n-1}(-z)^n}{k!(n-k)!}}$$
Shown here without Bell polynomials. Additionally, $k$’s upper bound can go up to $\infty$ with interchangeable sums in that case.. The radius of convergence uses Lambert W$(z)$: $|z|<\text W\left(\frac 1e\right),z\in\Bbb C$.
Switching the sums and limiting the index uses the lower regularized gamma function P$(a,z)$:
$$y=z-\frac12\lim_{c\to 0}\sum_{n=c,1+c,2+c,\dots}\frac{(-n z)^n\text P(2-n,-nz)}{e^{nz}nn!}$$
Shown here
