I'm trying to solve a heat transfer PDE in an annulus. This is the problem formulation: $$\frac{\partial u}{\alpha*\partial t} = \frac{1}{r} \frac{\partial u}{\partial r} + \frac{\partial^2 u}{\partial r^2}$$ The domain is: $ r_i < r < r_e $.
The annulus is subjected to a constant Temperature on $r_i$ and convection on $r_e$.
So, my B.C. are: $$ u(r_i,t) = T_0 $$ $$ \frac{\partial u(r_e,t)}{\partial r} + (h/k)*u(r_e,t) = (h/k)*T_r $$
To be consistent with the B.C. I chose this I.C.: $$ u(r,0) = \frac{(T_r - T_0)}{(k/h)+r_e-r_i} * r + T_0 - \frac{(T_r - T_0)*r_i}{(k/h)+r_e-r_i} $$
So I transformed my PDE into a new problem using: $$ u(r,t) = u^*(r,t) + w(r,t) $$ $$ w(r,t) = \frac{(T_r - T_0)}{(k/h)+r_e-r_i} * r + T_0 - \frac{(T_r - T_0)*r_i}{(k/h)+r_e-r_i} $$
The new problem is:
$$\frac{\partial u^*}{\alpha*\partial t} = \frac{1}{r} \frac{\partial u^*}{\partial r} + \frac{\partial^2 u^*}{\partial r^2} + \frac{(T_r - T_0)}{r*(k/h)+r_e-r_i} $$
Subject to:
$$ u^*(r_i,t) = 0 $$ $$ \frac{\partial u^*(r_e,t)}{\partial r} + (h/k)*u^*(r_e,t) = 0 $$ $$ u^*(r,0) = 0 $$
I hope everything up to now is correct. Then, using separation of variables with $ u^*(r,t) = R(r)*T(t) $ I obtained: $$ T(t) = C*e^{-\lambda^{2}\alpha t} $$ and a Bessel Equation: $$ \frac{\partial^{2} R(r)}{\partial r^{2}} + \frac{\partial R(r)}{r \partial r} + R(r) \lambda^{2} = - \frac{(T_r - T_0)}{r*(k/h)+r_e-r_i} $$.
However, before solving the Bessel equation with the forcing term, I found the eigenvalues of the problem by applying the B.C. and this is the expression that I obtained:
$$ - \lambda * J_1(r_e \lambda) Y_0(r_i \lambda) + \lambda * Y_1(r_e \lambda) J_0(r_i \lambda) + (h/k) * J_0(r_e \lambda) Y_0(r_i \lambda) - (h/k) * Y_0(r_e \lambda) J_0(r_i \lambda) = 0 $$
Now I now that we have infinite solution to the homogenuos Bessel equation $ R_\lambda(r) = A_\lambda*J_0(\lambda r) + B_\lambda*Y_0(\lambda r) $ since the $\lambda$ are infinite.
After this, I started looking for a way to solve the Bessel Equation with the forcing term $ - \frac{(T_r - T_0)}{r*(k/h)+r_e-r_i} $.
I tried to find a solution using the method of variation of parameters and i obtained a solution of this kind:
$$ u^*_{p, \lambda} = \frac{2}{\pi \lambda} * \frac{(T_r - T_0)}{(k/h)+r_e-r_i} * [ J_0(\lambda r) \int^{r_e}_{r_i} Y_0(\lambda r) dr - Y_0(\lambda r) \int^{r_e}_{r_i} J_0(\lambda r) dr]$$
(I used the fact that the Wronskian $ W(J_0(\lambda r),Y_0(\lambda r)) = \frac{\pi \lambda}{2r}$).
I'm concerned about the fact that I may have done some errors when deriving this solution. So,is this solution correct? When I try to integrate it, I obtain very different results from when I solve the very same ODE $ \frac{\partial^{2} R(r)}{\partial r^{2}} + \frac{\partial R(r)}{r \partial r} + R(r) \lambda^{2} = - \frac{(T_r - T_0)}{r*(k/h)+r_e-r_i} $ numerically with Matlab and I cannot understand where is the error; moreover, when using the solution $u^*_{p, \lambda}$ my B.C. are NOT satisfied.
It would be great if someone could check the derived solution up to this point and help me in solving the Bessel Equation with forcing term correctly. Thank you in advance.