Thank you topsquark for helping me.
The solution is very long and somehow complicated. I am not going to write every step, rather will be focusing on some important parts of it where the magic happens.
I am sure that you remember how to solve the Airy Equation using the power series [imath]\sum_{n=0}^{\infty}a_nx^n[/imath].
Solving the Bessel equation is similar but the solution will look like [imath]\sum_{n=0}^{\infty}a_nx^{n+m}[/imath] where [imath]m[/imath] is the indicial exponent, something related to Frobenius series if you recall that. It happens that we look for a solution including the indicial exponent in the power series due to regular singular points issues.
With this new form of the power series, we solve normally, taking the first and second derivatives, and substiting the results in the original Bessel differential equation.
After that, we get the recurrence relation
[imath]a_n = -\frac{1}{n(2\lambda+n)}a_{n-2}[/imath]
From that, we see that odd coefficients equal zero ([imath]a_1 = a_3 = a_5 = ........ = 0[/imath])
And we focus on even coefficients
[imath]a_2 = -\frac{1}{1 \cdot (\lambda + 1)} \cdot \frac{1}{2^2}a_0[/imath]
With a general formula this becomes,
[imath]a_{2n} = -\frac{1}{n \cdot (\lambda + n)} \cdot \frac{1}{2^2}a_{2n - 2}[/imath]
Don't ask me how, but this magically becomes,
[imath]a_{2n} = \frac{(-1)^n}{n!(\lambda + 1)(\lambda + 2).....(\lambda + n)}\cdot \frac{1}{2^{2n}}a_0[/imath]
With all this the solution becomes,
[imath]y(x) = a_0 x^{\lambda} + \sum_{n=1}^{\infty}a_{2n} x^{\lambda + 2n}[/imath]
After that, wihtout any further explanation, they say, the coefficient [imath]a_0[/imath] can have any nonzero value, and they chose,
[imath]a_0 = \frac{1}{2^{\lambda}\Gamma(1+\lambda)}[/imath]
And from here, the annoying Gamma function appeared. Wouldn't it be better to choose [imath]1[/imath] instead of Gamma function since [imath]1[/imath] is simple and nonzero value?
After that, the other coefficient magically becomes,
[imath]a_{2n} = \frac{(-1)^n}{2^{\lambda+2n} \ n! \ \Gamma(n + \lambda + 1)}[/imath]
[imath]y(x) = \frac{x^{\lambda}}{2^{\lambda}\Gamma(1+\lambda)} + \sum_{n=1}^{\infty}\frac{(-1)^nx^{\lambda + 2n}}{2^{\lambda+2n} \ n! \ \Gamma(n + \lambda + 1)} [/imath]
And this is in fact one of the independent solutions of the Bessel equation, namely, the the Bessel function of the first kind.
[imath]J_{\lambda}(x) = \frac{x^{\lambda}}{2^{\lambda}\Gamma(1+\lambda)} + \sum_{n=1}^{\infty}\frac{(-1)^nx^{\lambda + 2n}}{2^{\lambda+2n} \ n! \ \Gamma(n + \lambda + 1)} [/imath]
Okay, now I remember the gammas. You are talking about the coefficients. I had thought you were talking about as a function of x.
First, the recursion. Just work it through, n by n:
[imath]a_{2n} = \dfrac{-1}{4n(n+\lambda)} \cdot a_{2n - 2}[/imath]
From this we know that
[imath]a_{2n-2} = \dfrac{-1}{4(n-1)(n-1+\lambda)} \cdot a_{2n - 4}[/imath]
so putting that into the [imath]a_{2n}[/imath] equation
[imath]a_{2n} = \dfrac{-1}{4n(n+\lambda)} \left ( \dfrac{-1}{4(n-1)(n-1+\lambda} \right ) \cdot a_{2n - 4} [/imath]
and continuing:
[imath]a_{2n} = \dfrac{-1}{4n(n+\lambda)} \left ( \dfrac{-1}{4(n-1)(n-1+\lambda)} \right ) \left ( \dfrac{-1}{4(n-2)(n-2+\lambda)} \right )\cdot a_{2n - 6} [/imath]
Doing this n-1 times gives us
[imath]a_{2n} = \dfrac{-1}{4n(n+\lambda)} \left ( \dfrac{-1}{4(n-1)(n-1+\lambda)} \right ) \left ( \dfrac{-1}{4(n-2)(n-2+\lambda)} \right ) \dots \left ( \dfrac{-1}{4(1)(\lambda)} \right )\cdot a_0 [/imath]
or
[imath]a_{2n} = \dfrac{(-1)^n}{4^n n(n-1) \dots 1 \cdot (n + \lambda)(n + \lambda - 1) \dots \lambda}[/imath]
No magic to it. You should practice this. The solution of recursion relations isn't terribly well explained to Calculus students, but most of the time they aren't too bad to solve. You can usually spot the patterns in the derivations of the most common expansions pretty easily. (I'd advise you to learn how to solve at least 2nd order linear homogeneous recursions with constant coefficients. The theory is pretty similar to the corresponding ordinary differential equations. Inhomogeneous recursions can get nasty quickly, but there's a way to deal with those, too, though they are tougher than inhomogeneous differential equations.)
Now, [imath]n(n-1) \dots 1 = n![/imath], but here's where the problem comes in: nobody told us that [imath]\lambda[/imath] is a positive integer. The best we can do is say [imath](n + \lambda)(n + \lambda - 1) \dots \lambda = \Gamma(n + \lambda + 1)/\Gamma(\lambda + 1)[/imath]. There's no way around it unless we are lucky enough to get a whole number for [imath]\lambda[/imath] and, at least for a Physics problem, you can bet it won't be.
You get used to them after a while.
By the way, I'm not too keen on the rest of the derivation. We certain have the freedom to choose [imath]a_0[/imath] to be anything we want, and getting rid of one of the gamma functions is nice, but I think it's better to say we are going to absorb any constants into the arbitrary constants of the Bessel functions in the general solution form for y(x). It's simpler. (And that way we get rid of that [imath]2^{\lambda}[/imath]. I'm not sure why that's in there, either.)
-Dan
