Okay. An apology on two points:
1. We aren't doing the expansion into the Bessel series like I had thought. What we are really doing is
[imath]\displaystyle y(r,t) = C + \sum_{n=1}^{\infty} T_n(t) J_0( \sqrt{ \lambda_n } r)[/imath]
My only excuse is that I've never had to do that before and it didn't occur to me for this case. As it turns out, there should be a way to get it to work without the C term: the differential operator is the radial part of the Laplacian in cylindrical coordinates, which means we may think of y as a potential function. We can always add a constant to a potential and change nothing, so its presence or lack changes nothing in the end. But it does change how we define the expansion for f(r). It appears that your source did add it in and, frankly, it does make the solution look a tad nicer. The cost is that the [imath]f_n[/imath] aren't quite defined the way you might think.
2. Note that I dropped the [imath]a_n[/imath] in the above expansion. We really don't need it... the [imath]T_n(t)[/imath] will automatically take care of any extra constants. The [imath]a_n[/imath] just turns into an annoyance. Sorry about that.
[imath]\displaystyle \frac{dJ_n(x)}{dx} = \frac{1}{2}[J_{n-1}(x) - J_{n + 1}(x)][/imath]
I like the above identity that you used for the derivative of [imath]J_0[/imath].
I am just wondering why my identity isn't correct, although it gives the correct result?
The identity is
[imath]\dfrac{d}{dx} (x^n J_n(x)) = x^n J_{n-1}(x)[/imath]
I really don't know why you think you should be able to apply this to
[imath]\dfrac{d}{dx} \left ( x^n \dfrac{dJ_n}{dx} \right )[/imath]
As I said before, you need to apply the derivative rule first, simplify the [imath]J_{-1} = -J_1[/imath], then use the above identity. The reason you got the correct answer? Luck. You applied the identity incorrectly and it somehow came out right.
I'm going to add that extra C into your quote and make comments as needed.
I prefer first to apply the initial condition to find [imath]T(0),[/imath] so that I can solve the differential equation to find [imath]T(t)[/imath].
[imath]\displaystyle y(r,0) = C + \sum_{n = 1}^{\infty} a_n J_0( \sqrt{\lambda_n} r) T_n(0) = f(r)[/imath]
[imath]\displaystyle a_nT_n(0) = \frac{\int_{0}^{a} rJ_0(\sqrt{\lambda_n} r)f(r) \ dr}{||J_0(\sqrt{\lambda_n} r)||^2}[/imath]
Okay, let's stop a moment. The orthogonality condition is
[imath]\displaystyle \int_0^1 J_p ( \sqrt{\lambda_{pm}} r ) J_p ( \sqrt{\lambda_{pm}} r ) r \, dr = \dfrac{1}{2} \dfrac{1}{J_{p+1}^2(\sqrt{\lambda_{pn}})} \delta _{mn}[/imath]
So here we have
[imath]\displaystyle a_n T_n(0) = \dfrac{1}{2} \dfrac{1}{J_1^2(\sqrt{\lambda_n})} \int_0^1 (f(r) - C) J_0( \lambda_n r ) r \, dr[/imath]
Now all we need to do is equate this to our expansion:
[imath]\displaystyle C + \sum_{n=1}^{\infty} a_n T_n(t) J_0( \lambda_n r ) = f(r)[/imath]
[imath]\displaystyle \sum_{n=1}^{\infty} a_n T_n(t) J_0( \lambda_n r ) = f(r) - C[/imath]
and we may simply call the RHS
[imath]\displaystyle \sum_{n = 1}^{\infty} f_n J_0( \lambda_n r )[/imath]
which is what your source did. It's not
technically wrong to do this, but your source should have made clear that they were.
So equating the two series, we may immediately say that
[imath]a_n T_n(t) = f_n[/imath]
and we don't have to worry about the appearance of the C in the [imath]f_n[/imath] definition, and your error in the orthogonality condition gets erased.
Now you go on to apply the differential equation.
[imath]\displaystyle a_n\frac{dT_n}{dt} + a_n T_nk\lambda_n - q_n = 0[/imath]
[imath]\displaystyle T_n(t) = \frac{q_n}{a_nk\lambda_n} +c_1e^{-k\lambda_n t}[/imath]
where [imath]\displaystyle c_1 = T_n(0) - \frac{q_n}{a_nk\lambda_n}[/imath]
[imath]\displaystyle y(r,t) = C + \sum_{n = 1}^{\infty} a_nT_n(t)J_0( \sqrt{\lambda_n} r) [/imath]
[imath]\displaystyle y(r,t) = C + \sum_{n = 1}^{\infty} a_n\left(\frac{q_n}{a_nk\lambda_n} +c_1e^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} r) [/imath]
[imath]\displaystyle y(r,t) = C + \sum_{n = 1}^{\infty} \left(\frac{q_n}{k\lambda_n} +c_1a_ne^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} r) [/imath]
[imath]\displaystyle y(r,t) = C + \sum_{n = 1}^{\infty} \left(\frac{q_n}{k\lambda_n} +\left[T_n(0) - \frac{q_n}{a_nk\lambda_n}\right]a_ne^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} r) [/imath]
[imath]\displaystyle y(r,t) = C + \sum_{n = 1}^{\infty} \left(\frac{q_n}{k\lambda_n} +\left[a_nT_n(0) - \frac{q_n}{k\lambda_n}\right]e^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} r) [/imath]
Let [imath]f_n = a_nT_n(0)[/imath] <- Technically you are applying this, not "letting" it.
[imath]\displaystyle y(r,t) = C + \sum_{n = 1}^{\infty} \left(\frac{q_n}{k\lambda_n} +\left[f_n - \frac{q_n}{k\lambda_n}\right]e^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} r) [/imath]
We are almost arriving to the given solution. I will apply the boundary condition.
[imath]\displaystyle y(a,t) = C + \sum_{n = 1}^{\infty} \left(\frac{q_n}{k\lambda_n} +\left[f_n - \frac{q_n}{k\lambda_n}\right]e^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} a) = 0 \ \ [/imath] as [imath] \ \ J_0(\sqrt{\lambda_n} a) = 0[/imath]
Now this is taken care of.
but the boundary condition tells us that [imath]y(a,t) = T[/imath], so the only assumption that I can make is that I have to add a constant to our solution.
[imath]\displaystyle y(r,t) = C + \sum_{n = 1}^{\infty} \left(\frac{q_n}{k\lambda_n} +\left[f_n - \frac{q_n}{k\lambda_n}\right]e^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} r) [/imath]
[imath]\displaystyle y(a,t) = C + \sum_{n = 1}^{\infty} \left(\frac{q_n}{k\lambda_n} +\left[f_n - \frac{q_n}{k\lambda_n}\right]e^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} a) = C = T[/imath]
[imath]\displaystyle y(r,t) = T + \sum_{n = 1}^{\infty} \left(\frac{q_n}{k\lambda_n} +\left[f_n - \frac{q_n}{k\lambda_n}\right]e^{-k\lambda_n t}\right)J_0( \sqrt{\lambda_n} r) [/imath]
You did pretty good. The detail about the constant C wasn't your fault (your source apparently never mentioned it, and I gave you the wrong series expansion), and you did the rest reasonably well.
My suggestion is practice, practice, practice. Solving a differential equation by an orthogonal series expansion is actually a pretty basic technique, but some of the expansion functions themselves can get a bit... interesting. Bessel functions are one of the more interesting ones. I presume that you are still self-studying. I would suggest, after working some more problems with Bessel functions, you move on to Fourier series (if you haven't done them already), Legendre polynomials, Associated Legendre polynomials, Laguerre polynomials, and if you like you might at least glance at the hypergeometric functions. Those last are so general that practically every function you've seen up to now is actually a hypergeometric function, or directly related to one. I've never had to use them, myself (I rarely even use Bessel functions), so I wouldn't say that the hypergeometrics are critical, but I'm a Physicist and you are studying Mathematics: you might want to generalize a bit more than my advice.
-Dan
Addendum: The whole orthogonal expansion thing has a name: If you haven't seen it already, you might want to take a look at Sturm-Louiville theory. It's too advanced for your level, but it might be useful to take a look at it for future reference.