Solving 2-D steady state heat transfer in cylindrical coordinates I am trying to solve a 2-D steady state heat transfer equation in cylindrical coordinates
$$\frac{1}{r}\frac{\partial}{\partial r}\bigg(r\frac{\partial T}{\partial r}\bigg) +\frac{\partial^2 T}{\partial z^2}=0 ;\ 0\leq r \leq r_0; \ 0 \leq z \leq l$$
with BCs as follows:
$$T(r,0) = T_a$$
$$T(r,l)=T_h$$
$$\frac{\partial T}{\partial r}\bigg\rvert_{(0,z)}=\frac{\partial T}{\partial r}\bigg\rvert_{(r_0,z)}=0$$
What I have done so far:
Assume the solution of the form: $T(r,z)=R(r)Z(z)$
This lead to variable separation and a pair of ODEs
$$\frac{d}{d r}\bigg(r\frac{d R}{d r}\bigg) +\lambda^2rR=0;\  \frac{d^2 Z}{d z^2}-\lambda^2 Z=0$$
$\lambda^2$ being separation constant
Case 1 : (separation constant = 0)
This leads to $Z'' = 0$ and $rR'' + R' = 0$
The solution for $Z$ is $Z = A_1z+A_2$
The solution for $R$ is $R = C_1J_0(0) + C_2Y_0(0)$ ($J_0$ and $Y_0$ - Bessel functions of first and second kind respectively , order zero)
After applying BCs and the solution being bounded gives $$T = T_a + \bigg(\frac{T_h-T_a}{l}\Bigg)z$$
Case 2 (separation constant < 0)
This leads to $Z'' + \lambda^2 Z = 0 $ and $r^2R'' +rR' - \lambda^2 r^2 R =0$
The solution for $Z$ is $Z = A_1 \cos(\lambda z) + A_2 \sin(\lambda z)$
The solution for $R$ is $R = C_1 I_0(\lambda r) + C_2 K_0(\lambda r)$ ($I_0$ and $K_0$ - Modified Bessel functions of first and second kind respectively , order zero)
Applying BC at $r=0$ and realizing that the solution must be bounded here, $C_2 $ must vanish
Applying BC at $r=r_0$, restricts $\lambda$ to be zeros of $I_1(\lambda r_0)$
Hence the solution becomes $$T = \sum_{n=1}^{\infty} (A_n \cos(\lambda_n z)+B_n \sin(\lambda_n z))I_0(\lambda_n r)$$
Now, how do I find $A_n$ and $B_n$?
Case 3 : (separation constant > 0)
This leads to $Z'' - \lambda^2 Z = 0 $ and $r^2R'' +rR' + \lambda^2 r^2 R =0$
The solution for $Z$ is $Z = A_1 \cosh(\lambda z) + A_2 \sinh(\lambda z)$
The solution for $R$ is $R = C_1 J_0(\lambda r) + C_2 Y_0(\lambda r)$
Applying BC at $r=0$ and realizing that the solution must be bounded here, $C_2 $ must vanish
Applying BC at $r=r_0$, restricts $\lambda$ to be zeros of $J_1(\lambda r_0)$
Hence the solution becomes $$T = \sum_{n=1}^{\infty} (A_n \cosh(\lambda_n z)+B_n \sinh(\lambda_n z))J_0(\lambda_n r)$$
Now, how do I find $A_n$ and $B_n$?
 A: A lot of progress is made in the OP regarding the separable solutions of the ODE, in particular the rejection of solutions inconsistent with the insulating boundary conditions at $r=0,r_0$. (In particular, the solution set in the case of positive separation constant is on point.) However:


*

*When the separation constant is zero, the radial ODE $R''+rR'=0$ need not be solved in terms of Bessel functions; instead, one directly has $R(r)=1,r$ as the solution set. But the latter has $R'(0)\neq0$, so we reject it and keep only the $R(r)=1$ solution. Hence the relevant separable solution in this case is $T(r,z)=A_0+B_0 z$. (I'm not enforcing the boundary conditions in $z$ yet; see below).

*When the separation constant is negative, the radial solution set prior to enforcing the insulating BCs is indeed $R(r)=I_0(\lambda r),K_0(\lambda r)$. The OP also correctly notes that the second solution can be rejected on the basis that $\frac{d}{dr}I_0(\lambda r)$ vanishes at $r=0$ but $\frac{d}{dr}K_0(\lambda r)$ does not. However, $I_0(\lambda r)$ differs from $J_0(\lambda r)$ in that it does not oscillate for large $r$ but instead increases monotonically starting at $r=0$. So $\frac{d}{dr}I_0(\lambda r)=0$ has no solution besides $r=0$ and this solution set cannot satisfy both insulating BC's, and we must reject this case entirely.
With these corrections in mind, it should also be remembered that the final solution should be a linear combination of these separable solutions. Hence the final solution must be of the form
$$T(r,z)=A_0+B_0 z+\sum_{n=1}^\infty (A_n \cosh \lambda_n z+B_n \sinh \lambda_n z)J_0(\lambda_n r)$$
where $\lambda_n$ is the $n$th zero of $J_1(\lambda r_0)=0$.
At this point, we should enforce the boundary conditions at $z=0,l$. At $z=0$ we have $$T(r,0)=T_a= A_0+\sum_{n=1}^\infty A_n J_0(\lambda_n r).$$ But this will be true for all $0\leq r \leq r_0$ when $A_0=T_a$ and $A_{n\geq 1}=0$; moreover, this is the unique such solution. (I leave this as an exercise to the reader.) At $z=l$, one further obtains $B_0=\frac{T_h-T_a}{l}$ and $B_{n\geq 1}=0$, and this is again unique. From this we obtain the single unique result, which is the same one deduced by the OP in the case of zero separation constant. 
Note: From the result, we conclude that we didn't actually need any solutions other than those provided by the case of $\lambda=0$; they would, however, have been necessary had the boundary conditions at $z=0,l$ been functions of $r$ rather than constant.
