Is this expression for radial probability flux in Sakurai's Modern Quantum Mechanics wrong? The section on Schrodinger's equation for central potentials in Sakurai's Modern Quantum Mechanics (p. 208, 2nd edition) contains the following expression for the radial probability flux, as part of his argument for ruling out the (asymptotic) $r^{-l}$ solution for the radial part of the wavefunction, which I believe is wrong:
$$
j_r = \hat{\textbf{r}} \cdot \textbf{j} = \frac{\hbar}{m} R_{El}(r) \frac{dR_{El}(r)}{dr} \tag{1}
$$
Context:

Let us assume that the potential-energy function is not so singular so that  $\lim_{r\to0}r^2V(r) = 0$. Then, for small values of $r$, (3.7.9) [Radial part of Schrodinger equation as $r \to 0$] becomes
  $$
\frac{d^2u_{El}}{dr^2} = \frac{l(l+1)}{r^2}u_{El}(r)
$$
  [where $R_{El}(r)$ is the the radial part of the wavefunction and $u_{El}(r) = rR_{El}(r)$]
which has the general solution $u(r) = Ar^{l+1} + Br^l$
...
Consider the probability flux. This is a vector quantity whose radial component is 
  $$
j_r = \frac{\hbar}{m}Im(\psi^*\frac{\partial\psi}{\partial r}) = \frac{\hbar}{m}R_{El}(r)\frac{d}{dr}R_{El}(r)
$$

The textbook then goes on to substitute one by one, the asymptotic solutions $R = Ar^l$ and  $R = Br^{-(l+1)}$, and shows that for the second case, $lim_{r\to0}j_r \neq 0$.
Now, according to me, 
$$
\psi(r, \theta, \phi) = R_{El}(r)Y_l^m(\theta, \phi) \\
\frac{\partial \psi}{\partial r} = Y_l^m(\theta, \phi)\frac{dR_{El}(r)}{dr} \\
Im\left(\psi^*\frac{\partial \psi}{\partial r}\right) = \left|Y_l^m\right|^2Im\left(R^*\frac{dR}{dr}\right)
$$
For $R = Ar^{-(l+1)}$,
$$
R^*\frac{dR}{dr} = \left(A^* r^{-(l+1)} \right) \left(-(l+1)A r^{-(l+2)} \right) = -|A|^2(l+1)r^{-(2l+3)} \tag{2}
$$
which is purely real, thus for this case
$$
j_r = \frac{\hbar}{m} \left| Y_l^m \right|^2 Im\left(R^*\frac{dR}{dr}\right) = 0
$$
I looked up the errata (pdf) for Sakurai's book but the entry for page 208 only notes the absence of spherical harmonics from (1). Is there something wrong with my calculation (2)?
(Edited to add details of calculation)
 A: Consider l=0, so you don't think this has something to do with polar coordinates: perfect spherical symmetry, $\psi=R_E(r)\equiv \sqrt{\rho (r)} e^{iS(r)} $, for real ρ and S, in general. 
For bound states like atoms, of course, S vanishes, R is real, as you noticed, as E < 0 and these systems are stationary: they stay put,  $\dot{\rho}=0$, without leaking probability out, $j_r=0$. Atoms are stable.
But, for scattering states, of course, it does not: think of a spherical wave ($R\sim e^{irk}/r$, so $j_r=\hbar k/mr^2$), or a free spreading wave packet, below. 
The probability density current is thus 
$$
j_r= \frac{\hbar}{m}\operatorname{Im} \bar{R}\partial_r R= \frac{\hbar}{m} \rho \partial_r S .
$$
The continuity equation, $\dot{\rho}+ j_r(r)=0$ then implies that the probability P in a spherical volume V with surface area A, decreases as the current efflux from the surface shell,
$$
\dot{P}=\int dV ~\dot{\rho} = -\oint dA ~j_r. 
$$
So, for example, for the freely (zero is a spherical potential!) diffusing Gaussian wave packet, unnormalized, starting out with squared-width a at the origin of time (t =0),
$$
R=\left ( \frac{\sqrt{a}}{a+\frac{i\hbar t }{m}} 
   \right )^{3/2} ~\exp \left (-r^2 \frac{a-i\hbar t/m}{2(a^2+\hbar^2t^2/m^2)}  \right )~,$$ 
decidedly complex.
Evaluating the above current density,
$$
j_r= \frac{\hbar}{m} \rho \frac{r\hbar t }{m(a^2+\hbar^2 t^2/m^2 )}= \frac{r t \hbar^2}{m^2 a}\left ( \frac{1}{a+\frac{\hbar^2 t^2 }{m^2 a}} 
   \right )^{5/2} ~\exp \left (-r^2\frac{a}{a^2+\hbar^2t^2/m^2}  \right )~,
$$
you check the current and efflux attenuate with r (probability is conserved at the shell at infinity).  For $\hbar t/m \gg a$, a is replaced by $\hbar ^2t^2/m^2 a$, a width-squared growing to infinity, as the Gaussian collapses to a constant and localization is lost completely.
The "quantum flow velocity",
$$
\frac{j_r}{\rho}=\frac{r}{t}~\left ( \frac{1}{1+\frac{m^2 a^2}{\hbar^2 t^2}}\right ) 
$$
then collapses to r/t for large times.


*

*Edit on your amended question "where's my error in (2)?": Let us leave out the irrelevant spherical harmonics, etc, absorbing them as r-constants  in the term A the authors wish to banish. So, in my language, $\psi \propto R(r)=\sqrt{\rho}e^{iS(r)}\sim A/r^{l+1}$ near the origin. From my expression, $j_r\sim \frac{\hbar A A^*}{m r^{2l+2}}\partial_r S , $ so that the probability leakage from a shell near the origin is 
$$
4\pi r^2 j_r \propto \frac{\hbar A A^*}{m r^{2l}}\partial_r S \propto r^{-(2l+1)},$$
for a reasonable S vanishing at the origin, as you encounter in scattering theory. (Of course, if S is a constant, as in bound states, it is useless and absorbable, and the answer vanishes.) As a result of this efflux explosion, you have to reject this singular solution as unphysical. By contrast, for l=0, $S=kr$ is safe, amounting to a spherical wave, $R\sim \exp(irk) /r$. 

*A total aside, but while we are at it. For nonvanishing l and 
nonvanishing azimuthal q.n. m (not mass!), the bound state w.f. does have a complex behavior, $\exp(im\phi)$, leading to a nontrivial φ-component to the probability current! So, even though the probability never flows out of a spherical shell, there is this azimuthal probability flow like the jet stream, going round and round and round!  
