Special Differential Equation I ended up with a differential equation that looks like this:
$$\frac{d^2y}{dx^2} + \frac 1 x \frac{dy}{dx} - \frac{ay}{x^2} + \left(b  -\frac  c x - e x \right )y = 0.$$
I tried with Mathematica. But could not get the sensible answer. May you help me out how to solve it or give me some references that I can go over please? Thanks.
 A: Let $x=e^u$. I changed $e$ to $f$ in the equation to avoid confusions. Then, multiplying by $x^2$ gives $${x^2}\frac{{{d^2}y}}{{d{x^2}}} + x\frac{{dy}}{{dx}} - ay + \left( {b{x^2} - cx - f{x^3}} \right)y = 0$$
Now, if $x=e^u$, then $$\eqalign{
  & x\frac{{dy}}{{dx}} = \frac{{dy}}{{du}}  \cr 
  & {x^2}\frac{{dy}}{{dx}} = \frac{{{d^2}y}}{{d{u^2}}} - \frac{{dy}}{{du}} \cr} $$ so 
the equation  is
$$\frac{{{d^2}y}}{{d{u^2}}} - \frac{{dy}}{{du}} + \frac{{dy}}{{du}} - ay + \left( {b{e^{2u}} - c{e^u} - f{e^{3u}}} \right)y = 0$$
or $$\frac{{{d^2}y}}{{d{u^2}}} - ay + \left( {b{e^{2u}} - c{e^u} - f{e^{3u}}} \right)y = 0$$
$$\frac{{{d^2}y}}{{d{u^2}}} + \left( {b{e^{2u}} - c{e^u} - f{e^{3u}} - a} \right)y = 0$$
$$\frac{{{d^2}y}}{{d{u^2}}} + F\left( u \right)y = 0$$This is a $2^{\rm nd}$ degree DE. As Robert Israel points out, and shows in his answer, the solutions seem to be complicated. My inclination would be to aim for a series solution, finding the coefficients of $y$ recursively. 
A: Series solution around zero:
The point $x= 0$ is a regular singular. Taking the anzats
$$
y(x) = \sum_{n=0}^\infty q_n x^{n+s}
$$
we have
$$
\sum_{n=0}^\infty [(n+s)^2 - a]q_n x^{n+s-2} -\sum_{n=0}^\infty c q_n x^{n+s-1} + \sum_{n=0}^\infty b q_n x^{n+s} - \sum_{n=0}^\infty e q_n x^{n+s+1} =
$$
\begin{multline}
q_0(s^2 - a)x^{s-2} + \big[q_1\big((s+1)^2 - a\big) - c q_0\big] x^{s-1} + \\q_2 \big[\big((s+2)^2 - a\big) -c q_1 + b q_0\big]x^s + \sum_{n=0}^\infty \Big(...\Big)
\end{multline}
The indicial equation is $s^2-a=0$, hence $s = \pm\sqrt{a}$.
Case $s = \sqrt{a}$
In this case
\begin{align}
 q_1(2\sqrt{a} + 1) -c q_0 &=0,\\
 q_2(4\sqrt{a} + 2) -c q_1 + b q_0 &=0,
\end{align}
and the recurrence relation is
$$
 q_{m+3} = \frac{c q_{m+2} + b q_{m+1} - e q_m}{(m+3)(m + 3 + 2\sqrt{a})}.
$$
This, assuming all my algebra is right. 
What's important here is that
$$
 y(x) \sim x^\sqrt{a} z(x),
$$
which, combined with the form of $F(u)$ obtained by Peter Tamaroff, might help to propose a solution of the type
$$
 y(x) = x^\sqrt{a} \exp[\sqrt{F(u)} x] v(x),
$$
in a similiar way as it is done when solving the Hydrogen Atom or the Quantum Harmonic Oscillator as modern physics textbooks do (see Eisberg Fundamentals of Modern Physics, Hydrogen atom solution), that can lead to a relatively simple form for $v(x)$.

Edit 1
Taking the change of variables
$$
y(x) = \frac{z(x)}{\sqrt{x}},
$$
you end up with the equation
$$
z'' + \left\{\frac{1-4\alpha}{x^2} + \beta - \frac{\gamma}{x} - \epsilon x\right\}z = 0
$$
(where I've changed the constants to greek letters to avoid the $e$ confussion).
If $\epsilon = 0$ then, as Robert Israel points out, it reduces to the Whittaker Differential Equation (using the proper rescaling). Also, for $\alpha =1/4$ and $\gamma = 0$, you have the Airy Differential Equation.
Using the WKB approximation, 
$$
z(x) \sim A f^{-1/4} e^{\int f^{1/2} dx} + B f^{-1/4} e^{-\int f^{1/2} dx}
$$
where
$$
f(x) = \frac{1-4\alpha}{x^2} + \beta -\frac{\gamma}{x} -\epsilon x.
$$
you can try to find the asymptotic behavior of $z$, which I believe is
$$
z(x) \sim \frac{e^{-\frac{2}{3}(\epsilon^{1/3} x)^{3/2}}}{(\epsilon^{1/3} x)^{1/4}}
$$
for $x > 0$, and play with the differential equation resulting from taking
$$
z(x) = \frac{e^{-\frac{2}{3}(\epsilon^{1/3} x)^{3/2}}}{(\epsilon^{1/3} x)^{1/4}} v(x).
$$
A: I don't know if there  are closed form solutions in general.  In the case $e=0$, Maple finds a solution using Whittaker M and W functions:
$$y \left( x \right) =c_{{1}}
{{\rm \bf M}\left({\frac {ic}{2\sqrt {b}}},\,\sqrt {a},\,2\,i\sqrt {b}x\right)}
{\frac {1}{\sqrt {x}}}+c_{{2}}
{{\rm \bf W}\left({\frac {ic}{2\sqrt {b}}},\,\sqrt {a},\,2\,i\sqrt {b}x\right)}
{\frac {1}{\sqrt {x}}}
$$
Another interesting special case is $a=1/4$, $c=0$, where Maple's solution involves Airy functions:
$$
y \left( x \right) =c_{{1}}
{\text{Ai}\left(-{\frac {b-ex}{ \left( -e \right) ^{2/3}}}\right)}{\frac {1}{\sqrt {x}}}
+c_{{2}}{\text{Bi}\left(-{\frac {b-ex}{ \left( -e \right) ^{2/3}}}\right)}{\frac {1}{
\sqrt {x}}}
$$
EDIT: Note that the scaling $x \to k x$ preserves the form of the differential equation with $(a,b,c,e) \to (a,k^2b, kc,k^3e)$.  So if $e \ne 0$ we can assume WLOG that, say, $e=1$. 
As @Pragabhava noted, the indicial roots are $\pm \sqrt{a}$, so unless $\sqrt{a}$ is an integer there will be two fundamental solutions of the form 
$$\eqalign{y_1(x) &= x^{\sqrt{a}} \left(1 + \sum_{j=1}^\infty u_j x^j\right)\cr y_2(x) &= x^{-\sqrt{a}} \left(1 + \sum_{j=1}^\infty v_j x^j\right)}$$ with coefficients satisfying the recurrences $(2 \sqrt{a} j+j^2) u_j - c u_{j-1} + b u_{j-2} - u_{j-3} = 0$ (with $u_0 = 1$, $u_j = 0$ for $j < 0$) and 
$(-2 \sqrt{a} j+j^2) v_j - c v_{j-1} + b v_{j-2} - v_{j-3} = 0$ (with $v_0 = 1$, $v_j = 0$ for $j < 0$).  If $\sqrt{a}$ is an integer the second recurrence becomes singular at $j=2\sqrt{a}$, generally resulting in logarithmic terms.  I don't think there are closed-form solutions for the recurrences.
A: $\dfrac{d^2y}{dx^2}+\dfrac{1}{x}\dfrac{dy}{dx}-\dfrac{ay}{x^2}+\left(b-\dfrac{c}{x}-ex\right)y=0$
$\dfrac{d^2y}{dx^2}+\dfrac{1}{x}\dfrac{dy}{dx}-\left(ex-b+\dfrac{c}{x}+\dfrac{a}{x^2}\right)y=0$
Let $y=\dfrac{u}{\sqrt{x}}$ ,
Then $\dfrac{dy}{dx}=\dfrac{1}{\sqrt{x}}\dfrac{du}{dx}-\dfrac{u}{2x\sqrt{x}}$
$\dfrac{d^2y}{dx^2}=\dfrac{1}{\sqrt{x}}\dfrac{d^2u}{dx^2}-\dfrac{1}{2x\sqrt{x}}\dfrac{du}{dx}-\dfrac{1}{2x\sqrt{x}}\dfrac{du}{dx}+\dfrac{3u}{4x^2\sqrt{x}}=\dfrac{1}{\sqrt{x}}\dfrac{d^2u}{dx^2}-\dfrac{1}{x\sqrt{x}}\dfrac{du}{dx}+\dfrac{3u}{4x^2\sqrt{x}}$
$\therefore\dfrac{1}{\sqrt{x}}\dfrac{d^2u}{dx^2}-\dfrac{1}{x\sqrt{x}}\dfrac{du}{dx}+\dfrac{3u}{4x^2\sqrt{x}}+\dfrac{1}{x\sqrt{x}}\dfrac{du}{dx}-\dfrac{u}{2x^2\sqrt{x}}-\left(ex-b+\dfrac{c}{x}+\dfrac{a}{x^2}\right)\dfrac{u}{\sqrt{x}}=0$
$\dfrac{1}{\sqrt{x}}\dfrac{d^2u}{dx^2}-\left(ex-b+\dfrac{c}{x}+\dfrac{4a-1}{4x^2}\right)\dfrac{u}{\sqrt{x}}=0$
$\dfrac{d^2u}{dx^2}-\left(ex-b+\dfrac{c}{x}+\dfrac{4a-1}{4x^2}\right)u=0$
The above ODE is hypergeometric only when the special cases below:
$1$. $e=0$
$2$. $b=0$ and $c=0$
$3$. $c=0$ and $a=\dfrac{1}{4}$
Other than the above special cases the above ODE is not hypergeometric.
Unfortunately it also not belongs to any confluent forms of Heun’s equation.
Therefore to solve the above ODE generally is extremely difficult.
One of the main reason is that the coefficient of $u$ has too many terms or contains too high power terms. The similar situation also appear in Titchmarsh's ODE.
A: Consider a slightly more general ODE.
\begin{equation}
\frac{d^2 y(x)}{d x^2} + \frac{1}{x} \frac{d y(x)}{d x} + \left(-\frac{a}{x^2} + b - \frac{c}{x} - e x + e_1 x^2 \right) y(x)=0
\end{equation}
If we write:
\begin{equation}
y(x)=x^{\sqrt{a}}\cdot \exp\left( -\frac{\imath}{2 \sqrt{e_1}} x(-e+e_1 x)\right)  \cdot v(\frac{(-1)^{3/4}}{\sqrt{2} e_1^{1/4}} x)
\end{equation}
Then the function $v(x)$ satisfies the biconfluent Heun equation:
\begin{equation}
\frac{d^2 v(u)}{d u^2} -\left( \frac{\gamma}{u} + \delta + u\right)\frac{d v(u)}{d u} + \frac{\alpha u - q}{u} v(u) = 0
\end{equation}
where
\begin{eqnarray}
\delta&=& -1-2\sqrt{a}\\
\gamma&=& \frac{\left(\frac{1}{2}+\frac{i}{2}\right) e}{e_1^{3/4}}\\
\alpha&=& \frac{\imath \left(8 \imath \left(\sqrt{a}+1\right) e_1^{3/2}-4 b e_1+e^2\right)}{8 e_1^{3/2}}\\
q&=&\frac{\left(\frac{1}{4}+\frac{i}{4}\right) \left(2 \sqrt{a} e+2 i c \sqrt{e_1}+e\right)}{e_1^{3/4}}
\end{eqnarray}
 In[1304]:= Clear[y]; Clear[v]; Clear[m]; Clear[w]; Clear[f];
    m[x_] = Exp[-I/(2 Sqrt[e1]) x (-e + e1 x)] x^(Sqrt[ a]);
    y[x_] = m[x] w[x];
    myeqn = Collect[
       Simplify[(y''[x] + 
          1/x y'[x] + (-a/x^2 + b - c/x - e x + e1 x^2) y[x])], {w[x], 
        w''[x]}, Expand];
    myeqn1 = Collect[Simplify[myeqn/m[x]], {w[x], w'[x], w''[x], x^_}, 
       Simplify];

    T = (-1)^(3/4)/(Sqrt[2] e1^(1/4));
    f[x_] = T x;
    subst = {x :> f[x], 
       Derivative[1][w][x] :> 1/f'[x] Derivative[1][w][x], 
       Derivative[2][w][x] :> -f''[x]/(f'[x])^3 Derivative[1][w][x] + 
         1/(f'[x])^2 Derivative[2][w][x]};
    Collect[T^2 (myeqn1 /. subst /. w[f[x]] :> w[x]), {w[x], w'[x], 
      w''[x], x^_}, Simplify]


