Background
I am currently looking into the task of describing a transient temperature field $\theta(r,t)$ across the thickness $a \leq r \leq b$ of an infinitely long and hollow cylinder exposed to a sinusoidal temperature signal applied at the inner boundary $g(t)=\theta_0 \sin(2\pi \cdot f \cdot t)$.
The following thermal boundary conditions apply
$\theta(a,t)=g(t)$ $\quad$ for $t\ge0$, $\quad$ (1)
$\theta(b,t)=0$ $\quad$ for $t\ge0$, $\quad$ (2)
$\theta(r,0)=0$. $\quad$ (3)
The starting point for the derivation is the heat diffusion equation ($\alpha^*$ is the thermal diffusivity) in cylindrical coordinates
$\frac{\partial^2\theta}{\partial r^2} + \frac{1}{r}\frac{\partial\theta}{\partial r} = \frac{1}{\alpha^*}\frac{\partial\theta}{\partial t}$. $\quad$ (4)
The derivation of the final expression (which applies the Hankel transform followed by the inverse transform) for the temperature distribution can be found in several papers, for example [Shahani A.R. and Nabavi S.M. (2007) Applied Mathematical Modelling, Vol 31, p 1807-1818]. The final expression reads
$\theta(r,t)=-\alpha^*\pi \sum\limits_{n=1}^\infty \frac{\lambda^2_n J^2_0(\lambda_n b)}{J^2_0(\lambda_n a)-J^2_0(\lambda b)} \left[ Y_0(\lambda_n a)J_0(\lambda_n r) - J_0(\lambda_n a)Y_0(\lambda_n r) \right]\left[ e^{-\alpha^*\lambda_n^2 t} \int_0^t e^{\alpha^*\lambda_n^2\tau} g(\tau) d\tau \right]$ $\quad$ (5)
where $\lambda_n$ are the positive roots of the transcendental equation
$Y_0(\lambda_n a)J_0(\lambda_n b) - J_0(\lambda_n a)Y_0(\lambda_n b) = 0$ $\quad$ (6)
and $J_0(z)$ and $Y_0(z)$ are Bessel functions of order $0$ of the first and second kind, respectively.
Result
I have managed to implement the required code (as I understand it) in Matlab. The literature states that 100 roots is sufficient. When I plot the temperature field versus the thickness I get a strange temperature distribution for cases other than $g(t)=0$ as shown by the figure below (the blue and green curves are correct while the black and red curves are not).
My implementation does not seem to fulfill the boundary condition (1) which is also indicated by oscillations for the remainder of the corresponding (black and red) curves. My math skills are unfortunately rather limited and I do not seem to be able to see how the final expression (5) fulfills the boundary condition (1).
Question
Can it be shown that $Y_0(\lambda_n a)J_0(\lambda_n a) - J_0(\lambda_n a)Y_0(\lambda_n a) \ne 0$? Is this where my implementation goes wrong or is my mistake perhaps to be found elsewhere?