Non-self adjoint Sturm-Liouville problem I'm new to this site, but I felt the need to post when I recently came in to an ordinary differential equation/boundary value problem with this form:  
$(1)- \frac{d^2 y}{d x^2}  + \frac{m(m+1)} {x^2(1 + i \alpha e^{\beta x})} y = k^2 y$  
with the boundary conditions $y'(a) = 0$ and $ y'(b) = 0 $.  All the constant terms ${\alpha, \beta, m}$ are positive real numbers.  The differential operator induced, I'll call $L$, where  
$(2) L= - \frac{d^2 }{d x^2}  + \frac{m(m+1)} {x^2(1 + i \alpha e^{\beta x})} $ 
acts on the Hilbert space $L_2(a,b)$.  And  $0<a<x<b$, that is, $x$ defined on a positive real interval.  
The problem I'm experiencing stems from the fact that (1) is a non-self adjoint SL problem; where as self-adjoint problems are usually easier and the theory is well known.  I'm trying to prove existence and uniqueness of a solution to (1) and know as much as I can about its spectrum and (hopefully) find a rigorous or even a numerical solution (know the eigenvalues/eigenfunctions).  But, to be honest, I'm new to functional analysis and advanced differential equations need some advice.
So far, I've been able to write a formal solution to (1) using the Frobenius method, by rewriting (1) as 
$(3)-{x^2(1 + i \alpha e^{\beta x})} \frac{d^2 y}{d x^2}  + \left( m(m+1) - {x^2(1 + i \alpha e^{\beta x})}   k^2 \right) y = 0$  
and expanding the coefficients in (3) in a Taylor series at $x=0$.  This is not very useful though, the recursion formula is very messy and it's very difficult to tell where this solution converges.  Another option is to note that the coefficient of $y$ in (1) is infinitely differentiable and the derivatives can be written recursively.  That is, I can expand (1) in a Taylor series as well and use Frobenius, and get a different series to work with, which will likely not be much easier.  But that's again going to result in a messy recursion relation and (I presume) not be useful.
At this point, I'm a little lost and would really appreciate some professional advice.  Again, I'm looking for an existence/uniqueness proof, some theorems or analysis regarding the spectrum and a method for solution (numerical is okay) for the eigenvalues/functions other than what I've mentioned (admittedly, a tall order).  I should probably point out that this equation was the result of a separtion of variables of a partial differential boundary value problem in spherical coordinates.  The original PDE looked similar to the Helmholtz equation.  That's why the ODE above looks similar to the spherical version of Bessel's equation (e.g., set $\alpha=0$).
Regards.
 A: My first remark would be that this equation doesn't look like one whose solutions could be expressed with some standard special functions, so that way of expressing solutions is out of the question. However, for $\alpha = 0$ it looks like a version of the Bessel equation.
However, since your interval $[a,b]$ excludes $x=0$, all the coefficients are regular on it. Then the existence, uniqueness and regularity of solutions is standard: for any $c \in [a,b]$, the specifying the values of $y(c)$ and $y'(c)$ defines a unique solution that is smooth everywhere on $[a,b]$. Your coefficients are actually analytic in $x$, so the solution $y(x)$ will be analytic as well. Note also that your coefficients (except the leading one) also possibly depend analytically on the parameter $k^2$. In fact, those coefficients that depend on $k^2$ are entire functions of $k^2$ (no singularities for any $k^2\in \mathbb{C}$). Thus, if the values of $y(c)$ and $y'(c)$ are specified as independent of $k^2$, the solution $y(x;k^2)$ will actually be an entire analytic function of $k$ as well. (Coddington & Levinson, Theory of ordinary differential equations, McGraw-Hill 1955).
Now we come back to the boundary conditions. Let $y_a(x;k^2)$ be the solution satisfying $y(a)=1$ and $y'(a)=0$, and $y_b(x;k^2)$ be the solution satisfying $y(b)=1$ and $y'(b)=0$. Recall that for those values of $k^2$ that your differential operator $L-k^2$ is not invertible, are precisely your eigenvalues. If for some $k^2\in \mathbb{C}$ the two solutions $y_a(x;k^2)$ and $y_b(x;k^2)$ are independent then the operator $L-k^2$ is invertible and its inverse is an integral operator with the kernel (Green function)
$$
  (L-k^2)^{-1} = G(x,\xi;k^2) = \frac{y_a(x_<;k^2) y_b(x_>;k^2)}{W_{k^2}} ,
$$
where $x_< = \min(x,\xi)$, $x_> = \max(x,\xi)$ and $W_{k^2} = y'_a(x;k^2) y_b(x;k^2) - y_a(x;k^2) y'_b(x;k^2)$ is the Wronskian, which can be showed to be constant with respect to $x$. The formula may be off by a sign.
In spectral theory, the Green function (or rather the corresponding integral operator) would be referred to as the resolvent. The spectrum then corresponds precisely to the singularities of the resolvent in the $k^2$-complex plane. One can then go quite a long way in connecting the various properties of the above formula for $G(x,\xi;k^2)$ with the properties of the spectrum of $L$. The following monograph goes into depth in this direction, specifically concentrating on non-self adjoint boundary value problems: Mennicken & Möller, Non-Self-Adjoint Boundary Eigenvalue Problems, North-Holland 2003.
A: I made a little progress on this, and I am curious as to what people think.  The whole motivation for this approach is that it seems much easier to solve a self-adjoint problem than a non-self adjoint one.  There is so much more literature and development in the self-adjoint case, that I thought it may be worth while to try this:
Take the solution $y(x)$.  It's a complex function of a real variable.  We should be able to write it as $ y(x) = z(x) +i w(x)$  Where $z, w$ are purely real functions.  Also, rewrite the eigenvalue as $k^2= \Lambda + i \lambda$.  Now plug these back into the differential equation and separate the purely real and purely imaginary parts of the equation.  The result will be two coupled equations.  Like so.
$(4) -z'' +  \frac{m(m+1)} {x^2(1 + ( \alpha e^{\beta x})^2)} z +
 \frac{m(m+1) \alpha e^{\beta x} } {x^2(1 + ( \alpha e^{\beta x})^2)} w = \Lambda z - \lambda w $
$(5) -w'' +  \frac{m(m+1)} {x^2(1 + ( \alpha e^{\beta x})^2)} w -
 \frac{m(m+1) \alpha e^{\beta x} } {x^2(1 + ( \alpha e^{\beta x})^2)} z = \Lambda w + \lambda z 
$
Now, looking at (4) and (5), two things seem possible.  One, for example, is to solve (5) algebraically for $z$ and plug that into (4).  In which case the result will be a fourth order self-adjoint problem.  
The other possibility is to view (4) as an inhomogeneous equation for $z$, where $\Lambda$ plays the part of the eigenvalue and the terms with $w$ act as the inhomogeneous "forcing" term.  Similar remarks apply to equation (5).  In this view the "homogeneous" form of (4) is a self-adjoint second order Sturm-Liouville problem and so is (5).  Therefore, $w,z$ can be exapanded in terms ortho-normal eigenfunctions.
It seems then, that one would have to solve for the $\Lambda$ eigensystem, expand $w,v$ in terms of those eigenfunctions, then use this somehow to get the $\lambda$.  I'm still a bit sketchy on a few things, and appreciate any comment on this.
A: There is theory for "accretive"  (= minus "dissipative") operators in Hilbert space, which have numerical range in a half-plane (also "sectorial" operators). See Kato's Perturbation Theory
for Linear Operators Section V.3.10, p. 278 and VI.1.5 p. 318. I think you should have stability of eigenvalues and eigenfunctions for small $\alpha$ and $\beta$. (Also, VII Analytic perturbation theory p.364. I don't know enough to say if this is applicable.)
