The Hamiltonian for the spin-spin interaction is:
$\Delta H_{SS} = \frac{\gamma_p e^2}{m m_p c^2 r^3} \Big( \frac{1}{r^3} \big(3(\vec{s}_p \cdot \hat{r})(\vec{s}_e \cdot \hat{r})-(\vec{s}_p \cdot \vec{s}_p) \big)+\frac{8 \pi}{3} (\vec{s}_p \cdot \vec{s}_p) \delta^{(3)}( \vec{r} ) \Big) $
For the cases where $l \neq 0$ the term with the delta function cancels, and the wavefunction is proportional to $r^l$ for small $r$ values. Thus, when $l>0$ we get that $\psi(0)=0$, and then the correction to the energy will be:
$\Delta E_{hf} = \frac{\gamma_p e^2}{m m_p c^2} \langle \frac{1}{r^3} \big( ( \vec{l} \cdot \vec{s}_p )+3(\vec{s}_p \cdot \hat{r})(\vec{s}_e \cdot \hat{r})-(\vec{s}_p \cdot \vec{s}_p) \big) \rangle $
This expectation value was calculated by Bethe and Salpeter, and the result is:
$\Delta E_{hf} = \frac{m}{m_p} \alpha^4 mc^2 \frac{\gamma_p}{2 n^3} \Big( \frac{ f(f+1)-j(j+1)-\frac{3}{4}}{j(j+1)(l+\frac{1}{2})} \Big)$ ,
This result coincide with the $l=0$ case, since then $j=\frac{1}{2}$ and the proton has spin 1/2 so $f=j \pm \frac{1}{2}$ and the expression above can then be simplified to the result you already have for the $l=0$ case...
(Next time try older QM books like this one http://adsabs.harvard.edu/abs/1957qmot.book.....B )