I figured out that the problem is that in your code you are using besselj in MATLAB which is ordinary bessel function of first kind (J_n). You should have used spherical bessel functions j_n instead in this expansion.
The spherical Bessel functions j_n are given in terms of ordinary Bessel...