Eigenvalues and eigenfunctions of the exponential potential $ V(x)=\exp(|x|) $ For $a$ being positive what are the quantisation conditions for an exponential potential?
$$ - \frac{d^{2}}{dx^{2}}y(x)+ ae^{|x|}y(x)=E_{n}y(x) $$
with boundary conditions  $$ y(0)=0=y(\infty) $$
I believe that the energies $ E_{n} $ will be positive and real 
I have read a similar paper:

P. Amore, F. M. Fernández. Accurate calculation of the complex eigenvalues of the Schrödinger equation with an exponential potential. Phys. Lett. A 372 (2008), pp. 3149–3152, arXiv:0712.3375  [math-ph].

However, I get this strange quantisation condition        
$$ J_{2i\sqrt{E_{n}}}(\sqrt{-a})=0 $$
How can I handle this in the case $ a >0 $?
 A: The paper you quote covers a similar case, which was solved previously by S.T. Ma (Phys. Rev. 69 no. 11-12 (1946), p. 668), but deals with the scattering problem on the tail of the exponential - hence the complex energies. What follows is partly inspired by that paper but is quite distinct from it. The tricky part is not getting scared by the Bessel functions, but that's why we have the theory of special functions.
For one, the exponential potential $e^{|x|}$ you ask for is an even function, which means that the corresponding eigenfunctions on $(-\infty, \infty)$ will be either even or odd. Therefore, they can be treated as eigenfunctions of the simpler potential $e^x$ on $(0,\infty)$, with boundary condition $\psi'(0)=0$ or $\psi(0)=0$, respectively. Since you ask about the latter condition, there is no point in keeping the absolute value.
The problem, then, is the eigenvalue problem
$$-\frac{d^2}{dx^2}\psi+A e^x \psi=E\psi\text{ under }\psi(0)=0=\psi(\infty).\tag{problem}$$
(A word on dimensions: To get the equation down to this form, we've had to set $\hbar, m$ and the length scale of the exponential to 1, by taking appropriate units of time, mass and length respectively. This means that there is no more dimensional freedom, and the hamiltonian has one free parameter, $A$, which will affect not only the scale of the spectrum (which you might expect as $A$ is a scaling on the potential) but also its structure.)
Amore et al. treat this as a boundary-value problem in $\mathbb C$ and using a change to a complex variable. This complicates the issue more than is really necessary and for simplicity I will use only real variables, though this comes at the cost of dealing with modified Bessel functions instead of standard ones. The initial step is to change variable to $z=2\sqrt{A}e^{x/2}$, so that $Ae^x=z^2/4$ and derivatives transform as
$$
\frac {\partial }{\partial x}=\frac {\partial z}{\partial x}\frac {\partial }{\partial z}=\frac {z }{2}\frac {\partial }{\partial z}
\text{ so }
\frac {\partial^2 }{\partial x^2}
=\frac14\left(
z^2\frac {\partial^2 }{\partial z^2}+z\frac {\partial }{\partial z}
\right).
$$
The final equation is thus
$$
\left[
z^2\frac {\partial^2 }{\partial z^2}+z\frac {\partial }{\partial z}-(z^2+\nu^2)
\right]\psi=0
\tag{equation}
$$
where $\nu=i\sqrt{4E}$. (Yes. Some complexness is inevitable. No fear, it will eventually not matter.)
This equation is Bessel's equation in modified form with index $\nu$. This is exactly the same as Bessel's equation for more normal situations; the index is complex but that is all. Two linearly independent solutions are the modified Bessel functions of the first and second kind, $I_\nu(z)$ and $K_\nu(z)$, so the general solution of $(\text{equation})$ looks like
$$
\psi(z)=aI_\nu(z)+bK_{\nu}(z).
$$
We then only need to impose the boundary conditions $\psi|_{z\rightarrow \infty}=\psi|_{z=2\sqrt{A}}=0$:


*

*The condition at infinity requires that we set the coefficient of $I_\nu$ to zero, since the First Kind function always explodes. We could have done this from the start: $K_\nu$ is, by definition, the exponentially decaying solution, while $I_\nu$ grows exponentially.

*The condition at $x=0$ then simply requires that $K_\nu(2\sqrt{A})=0$. In terms of energies, then,
$$
\boxed{K_{2i\sqrt E}(2\sqrt{A})=0,}
$$
and this is your quantization condition.
As it happens, $K_\nu(z)$ is real for real $z$ and purely imaginary $\nu$. One way to prove this is via this integral representation:
$$
K_{\nu}(x)
=\sec(
{\nu\pi}/{2})\int_{0}^{\infty}\cos(x
\sinh\nolimits t)\cosh(\nu t)dt,
$$
which is the analogue of Bessel's First Integral for $K_\nu$; alternatively, it also follows from the conjugation symmetry
$$
K_{\nu^*}(z^*)=K_\nu(z)^*,
$$
which follows from the integral representation but is a more fundamental property of the function.
Since $K_\nu$ is real here, for whatever reason, we can ask for its zeros. As with all Bessel zeros there is no chance of an elementary formula for them, but they can be found quite easily using numerical methods (for properties of the zeros, see this DLMF reference). For a taster, here are some graphs, in log-linear scale (so zeros show up as downward, log-like peaks), of $K_{2i\sqrt{E}}(2\sqrt{A})$ as a function of $E$, for a few different values of $A$.

While there isn't all that much to say about the energies from this, it is clear that there are a countable infinity of them, that they are bigger than $A$, and that their spacing increases with increasing $A$ and $n$ (why?) - but that's really all you'd really want to know!

Just for completeness: the eigenfunctions themselves, then, are of the form
$$\psi_n(x)=C_nK_{2i\sqrt{E_n}}\left(
2\sqrt{A}e^{x/2}
\right).$$
It is interesting to note that the dependence in $n$ comes through the index instead of a coefficient before $x$. This is partly to ensure the very strict decay $\psi\sim e^{-\exp(x/2)}$, which is required by the very hard exponential wall of the potential. For some information on how these Bessel functions behave, try the Functions of Imaginary Order subsection in the DLMF; particularly important results are asymptotics on $K_{i|\nu|}$ at large $x$ and for the oscillatory region. The latter is 
$$
{K}_{{i\nu}}\left(z\right)=-\left(\frac{\pi}{%
\nu\mathop{\sinh}\left(\pi\nu\right)}\right)^{{\frac{1}{2}}}%
\mathop{\sin}\left(\nu\mathop{\ln}\left(\tfrac{1}{2}%
z\right)-\gamma_{\nu}\right)+\mathop{O}\left(z^{2}\right),
$$
so the asymptotic for the wavefunction is of the form $\psi(x)\sim\sin\left(\sqrt{E_n}x\right)$, as it should be. (Note, though, that this holds little physics beyond the standard: the information on the potential's variation is encoded in the change of the instantaneous frequency as in e.g. these formulas, and would require beefier maths.)
