I'm trying to obtain the analytical solution of a Fokker-Planck PDE, which the solution is a probability density function, and then use this to find the mean of some quantity in the paper. The paper has a solution which they say can be found via Mehler-Fock transform. Their solution which I am trying to obtain reads
$$P_{\varepsilon}(L,u) = \frac{e^{-\varepsilon^2sL/4}}{2\sqrt{2\pi}(\varepsilon^2sL)^{\frac32}}\int_{u}^{\infty}\frac{xe^{-x^2/(\varepsilon^2sL)}}{\sqrt{\cosh{(x)}-\cosh{(u)}}}dx.$$
I obtain a different solution to this one.
$\textbf{My attempt:}$ The equation in the paper reads \begin{align} \frac{\partial P_\varepsilon}{\partial L} = \varepsilon^2s\frac{\partial}{\partial u}\bigg[(u^2-1)\frac{\partial P_\varepsilon}{\partial u}\bigg], \quad P_\varepsilon(L=0,u) = \delta(u-1). \end{align} This is a $1d$ Fokker-Planck equation, where $P_\varepsilon(L,u)$ is a probability density for some diffusion Markov process. To solve this equation, note that the right hand side is simply a Legendre differential equation. Denote the legendre function of the first kind via \begin{align} \frac{d}{du}(u^2-1)\frac{d}{du}P_{-\frac12+i\mu}(u) = -\bigg(\mu^2+\frac14\bigg)P_{-\frac12+i\mu}(u), \end{align} which has an integral representation \begin{align} P_{-\frac12+i\mu}(u) = \frac{\sqrt{2}}{\pi}\cosh{(\pi\mu)}\int_{0}^{\infty}\frac{\cos{(\mu\tau)}}{\sqrt{\cosh{(\tau)}+u}}d\tau. \end{align} Now use the Mehler-Fock transform. The Mehler-Fock transform of an integrble function $f$ defined on $[1,\infty)$ is the function $\check f$ defined on $[0,\infty)$ where \begin{align} \hat{f}(\mu) = \int_{1}^{\infty}f(u) P_{-\frac12+i\mu}(u)\,du, \end{align} with inverse transform \begin{align} f(u) = \int_{0}^{\infty}\check{f}(\mu)\mu\tanh{(\mu\pi)}P_{-\frac12+i\mu}(u)\,d\mu. \end{align} Applying the Mehler-Fock transform to the PDF $P_\varepsilon(L,u)$ gives \begin{align} \check p(L,\mu) = \int_{1}^{\infty}p(L,u)P_{-\frac12+i\mu}(u)du. \end{align} Taking a partial derivative in $L$ gives \begin{align} \frac{\partial\check p}{\partial L}(L,\mu) = \varepsilon^2s\int_{1}^{\infty}\frac{\partial}{\partial u}\bigg[(u^2-1)\frac{\partial p}{\partial u}(L,u)\bigg]P_{-\frac12+i\mu}(u)du. \end{align} Integrating twice more by parts gives \begin{align} \frac{\partial\check p}{\partial L}(L,\mu) = \varepsilon^2s\int_{1}^{\infty}p(L,u)\frac{\partial}{\partial u}\bigg[(u^2-1)\frac{\partial P_{-\frac12+i\mu}}{\partial u}(u)\bigg]du. \end{align} Using the ODE $(2)$ which satisfies the Legendre function, the Mehler-Fock transform satisfies the ODE \begin{align} \frac{\partial \check p}{\partial L}(L,\mu) = -\varepsilon^2s\bigg(\mu^2+\frac14\bigg)\check p(L,\mu), \quad \check p(L=0,\mu) = 1. \end{align} Then \begin{align} \check p(L,\mu) = \exp{\bigg(-\bigg(\mu^2+\frac14\bigg)L\varepsilon^2s\bigg)}. \end{align} Hence, the solution to $(1)$ is \begin{align}\nonumber P_\varepsilon(L,u) &= \int_{0}^{\infty}\mu\tanh{(\mu\pi)}P_{-\frac12+i\mu}(u)\exp{\bigg(-\bigg(\mu^2+\frac14\bigg)L\varepsilon^2s\bigg)}d\mu \\ &= \int_{0}^{\infty}\mu\tanh{(\mu\pi)}P_{-\frac12+i\mu}(u) = \frac{\sqrt{2}}{\pi}\cosh{(\pi\mu)}\int_{0}^{\infty}\frac{\cos{(\mu\tau)}}{\sqrt{\cosh{(\tau)}+u}} \\ &\times\exp{\bigg(-\bigg(\mu^2+\frac14\bigg)L\varepsilon^2s\bigg)}d\tau d\mu. \end{align}
This paper claims to have a solution
$$P_{\varepsilon}(L,u) = \frac{e^{-\varepsilon^2sL/4}}{2\sqrt{2\pi}(\varepsilon^2sL)^{\frac32}}\int_{u}^{\infty}\frac{xe^{-x^2/(\varepsilon^2sL)}}{\sqrt{\cosh{(x)}-\cosh{(u)}}}dx.$$ Here I have computed my solution vs theirs, for varying parameter $L\in[1,10]$ and fixed $\varepsilon^2s$
Clearly they do not agree.
$\textbf{My second question}$ is how they use this solution to find a mean value via
\begin{align}\mathbb{E(R)} &= \int_{1}^{\infty}\bigg(\frac{u-1}{u+1}\bigg)P_{\varepsilon}(L,u)du \\ &= \frac{e^{(-\varepsilon^2sL/4)}}{2\sqrt{2\pi}(\varepsilon^2Ls)^{\frac32}}\int_{1}^{\infty}\bigg(\frac{u-1}{u+1}\bigg)\int_{u}^{\infty}\frac{xe^{-x^2/(4\varepsilon^2sL)}}{\sqrt{\cosh{(x)}-\cosh{(u)}}}dxdu \end{align}
$\textbf{AND then reduce this expression to read}$
$$\mathbb{E(R)}=1-\frac{4}{\sqrt{\pi}}e^{-\varepsilon^2sL}\int_{0}^{\infty}\frac{x^2e^{-x^2}}{\cosh{(\sqrt{\varepsilon^2sL}x)}}dx,$$ since it's easy to see asymptotically that they behave differently. I'm looking for help in either of these questions. Thanks for the help in advance.