ODE with singular coefficients I started with an ODE (first ODE) :
$-(1-x^2)y''(x) +x y'(x) -  q(x) y(x) = \lambda y(x).$ 
Then I got a more sophisticated differential equation ( second one) 
and is given by
$$-(1-x^2)y''(x) +x y'(x) -  q(x)y(x) + \frac{1}{(1-x^2)}y(x) = \lambda y(x).$$ 
What I want is to understand the following: 
Given a solution of equation (1), how can I generalize this one to a solution of equation (2). 
 A: With convenient change of variable, the ODE can be transformed to another ODE on a simplified form :  $\frac{y’’(t)}{y(t)}=P_2(\cosh(t))$ where $P_2$ is a polynomial of the second degree.
It is known that some particular ODE’s of this form can be solved in terms of MathieuS and MathieuC functions: http://mathworld.wolfram.com/MathieuDifferentialEquation.html
For example, this is the case if one of the coefficient is null and the two others not null : in the case $\gamma=0$ or in the case $\alpha=0$. But not in the case $\lambda=0$. Other cases are solvable in case of particular relationship between the coefficients, in terms of Mathieu funcions or in terms of other special functions of lower level for particular values of the coefficients of the polynomial. 
In the general case, i.e. $\alpha, \gamma, \lambda$ not null and not related one to another, a closed form for the solutions would certainly require special functions of higher level, which are not referenced as standard functions. So, I think that the ODE on its general form cannot be analytically solved in that sense.

The second EDO with one more term becomes on the reduced form :

A fortiori, in the general case, the solutions cannot be expressed with a finite number of standard functions. But closed forms of solutions could exist in some cases of particular values of the coefficients $\alpha, \gamma, \lambda, \beta$ 
Adding a term of different nature to an ODE generally changes a lot the form of the solution. A simplified example is shown below : The solution of the ODE which was expressed in terms of elementary functions becomes much more complicated because requiring special functions (associated Legendre) :

Another example which involves even more complicated functions :

A solution on the exponential form exists in case of particular relationships between the parameters $\gamma, \alpha, \lambda$ :

With the additional term, an approximative solution on exponential form can also be found  in case of particular relationships between the parameters $\gamma, \alpha, \lambda, \beta$ :
 
But this is of low interest for practical application because it cannot be generalized to ODEs which parameters are not related on the same way.
All this draw to think that the analytic method is not convenient to tackel this problem in case of a so complicated ODE. Other approaches might be more convenient, especially systematic studies with numerical methods.
A: You can transform your final equation to Sturm-Liouville form by the substitution
$$
                y = (1-x^{2})^{1/4}f.
$$
I believe that the new equation has the Sturm-Liouville form
$$
         -\frac{d}{dx}(1-x^{2})\frac{df}{dx}+\left[\frac{u^{2}}{1-x^{2}}+\frac{1}{4}-\alpha x-\gamma x^{2}\right]f= \lambda f.
$$
This form makes the equation look a little more classical. And it seems to suggest you'll have an orthogonality for the original solutions in the form
$$
           \int_{-1}^{1}\frac{1}{\sqrt{1-x^{2}}}y_{1}y_{2}\,dx = 0.
$$
Of course you'll need to find the right endpoint conditions.
Classical Case: If you start with the equation
$$
        (1-x^{2})y''-2xy'+n(n+1)y = 0,
$$
and you differentiate all terms $m$ times, you get
$$
    (1-x^{2})y^{(m+2)}+{{m}\choose{1}}(-2x)y^{(m+1)}+{{m}\choose{2}}(-2)y^{(m)}\\
    -2xy^{(m+1)}+{{m}\choose{1}}(-2)y^{(m)}\\+n(n+1)y^{(m)}=0.
$$
That gives you a different form of the Associated Legendre Equation
$$
    (1-x^{2})y^{(m+2)}-2(m+1)xy^{(m+1)}+[n(n+1)-m(m-1)-2m]y^{(m)}=0. (?)
$$
This is a related second-order differentiation equation for $y^{(m)}$ which, after a change of variable $y=(1-x^{2})^{-m/2}f$ leads to the standard form of the Associated Legendre equation.
Even if you try clearing denominators in your posted DiffEQ and trying the same trick, the order of the derived equation will be too high. The classical trick works only because the coefficients are polynomial with orders that match the order of the derivative for that term, including the lowest order term which is classically constant.
You need some other type of ladder operator, and I'm not good at constructing those devices. BTW: I think it was Schrodinger who came up with ladder operators. I ran across one of his early publications on that. Amazing genius.
A: This isn't an answer but it's too long for a comment. Let $t=x^2$ and note that 
$$ x \frac{d}{dx}=x \frac{dt}{dx}\frac{d}{dt}=x\cdot 2x\frac{d}{dx}=2t\frac{d}{dt},\\
   x^2\frac{d^2}{dx^2} = x \frac{d}{dx}\left(x \frac{d}{dx}-1\right)=2 t \frac{d}{dt}\left(2t\frac{d}{dt}-1\right)=4t^2 \frac{d^2}{dt^2}+2t \frac{d}{dt}.$$
If we then multiply your differential equation by $x^2$, we can thus use these identities to rewrite as $$-4t^2(1-t)y''(t)+2t(2t-1)y'(t) \\- 4 (2t-1) \left( \alpha  + \gamma (2t-1) \right) y(t)+\frac{4u^2-1}{4t(1-t)}y(t) = \lambda \, t \, y(t)$$
This may be easier to approach, though it's still quite busy.
