Expectation value of $p^2 (1/r) + (1/r) p^2$ I'm trying to derive $\langle kl | \lbrace 1/r , p^2 \rbrace | nl \rangle$, where the states satisfy the equation of motion (I omit factors of $1/2m$ etc.):
$$(p^2 + V)| n,l \rangle = E_n | n,l \rangle$$
and $\lbrace A,B \rbrace = AB + BA$ is the anticommutator.
I have two solutions at hand that differ. In general I can write:
$$\langle k,l | \lbrace 1/r , p^2 \rbrace | n,l \rangle = \langle k,l | 1/r~p^2 + p^2~1/r | n,l \rangle$$
and I can use that $p^2 \propto -\nabla^2$ and $\nabla^2 1/r \propto -\delta^{(3)}(r)$.
Now I want to let $p^2$ act first on the respective states and let $1/r$ act later on, such that I can use the e.o.m. Thus I have:
$$\langle k,l | 1/r~p^2 + p^2~1/r | n,l \rangle = \langle k,l | 1/r~(E_n - V) + (E_k - V)~1/r | n,l \rangle$$
In this solution I let the first summand act on the ket and I let the second summand act on the bra. However, when I do it the other way around, using the product rule $\nabla^2 (1/r~\psi) = (\nabla^2 1/r) \psi + 1/r (\nabla^2 \psi)$ and use the e.o.m for the second term I obtain:
\begin{eqnarray}
\langle k,l | 1/r~p^2 + p^2~1/r | n,l \rangle &=& \langle k,l | \delta^{(3)}(r) + (E_k - V) + (E_n - V) + \delta^{(3)}(r) | n,l \rangle\\ &=& \langle k,l | 1/r~(E_n - V) + (E_k - V)~1/r | nl \rangle + 2 \langle kl | \delta^{(3)}(r) | n,l \rangle
\end{eqnarray}
These two solutions are clearly different if and only if both of the wave functions do not vanish at the origin. This is to say that these two solutions give different results for s-waves.
Am I missing something essential here? Any input would be appreciated!
 A: A paradigm is worth a hundred puffs of bloviation.
Exploit your spherical symmetry, $\hat{r}\cdot \vec{\nabla}=\partial_r$, and $\vec{\nabla} f(r) =\hat{r} \partial_r f(r)$; and, since the challenge you are concerned with is the singularity at the origin, r =0, we might as well drop all l>0, which are softer than s-waves, as you noted. 
In your $\hbar=1$, $m=2$ units, consider, for simplicity, the s-waves of the Hydrogen hamiltonian, sufficient to illustrate the point,
$$
H=p^2-2/r= -\frac{1}{r^2} \partial_r (r^2\partial_r) -\frac{2}{r}~.\tag1
$$
The Hilbert space integrations will be over $4\pi \int^\infty_0 dr~ r^2 ~~$ .
Thus,
$$
\nabla^2 \frac{1}{r}=-4\pi \delta^{(3)} (\vec{r})~,\tag2
$$
which, in our radial context, amounts to just 
$$
-p^2 \frac{1}{r}=\frac{1}{r^2} \partial \left(r^2\partial \frac{1}{r}\right )=- \frac{\delta(r)}{r^2}~.\tag{2'}
$$
The ground eigenstate (n =1) is then
$$
\psi_0= \frac{e^{-r}}{\sqrt{\pi/2}}   \qquad \Longrightarrow  E_0=-1, \tag3
$$
while the first excited state (n =2) is 
$$
\psi_1=  \frac{e^{-r/2}}{4\sqrt{2\pi}}  (2-r)   \qquad \Longrightarrow  E_1=-1/4, \tag4
$$
so that $\langle \psi_1| \psi_0\rangle=0$, since $\int_0^\infty \! dx~ e^{-x} x^n =n! ~~$ .

It follows that 
$$
\frac{1}{r}~ p^2 ~\psi_0= \left (\frac{2}{r^2}-\frac{1}{r}\right )\psi_0 \tag5 \\
 p^2 \left (\frac{1}{r} \psi_0 \right) = \left(\frac{\delta(r)}{r^2}-\frac{1}{r}\right )\psi_0 . 
$$
(The second eqn is the screened Poisson eqn -- "massive photon".) Hence 
$$
[p^2,1/r]\psi_0 = \left ( \frac{\delta(r)}{r^2}-\frac{2}{r^2}  \right) \psi_0. \tag{5'}
$$
(In general, $[p^2,1/r]=\frac{\delta(r)}{r^2} + \frac{2}{r^2} \partial_r ~$ .)
Now note that 
$$
\langle \psi_0|[p^2,1/r]| \psi_0\rangle=0= \langle \psi_1|[p^2,1/r]| \psi_1\rangle .\tag6
$$
Consequently, as per your first evaluation path using hermiticity, 
$$
\langle \psi_0|\{p^2,1/r\}| \psi_0\rangle=
\langle \psi_0|(H+2/r)  \frac{1}{r} +  \frac{1}{r} (H+2/r)  | \psi_0\rangle=
\langle \psi_0| \frac{2}{r} \left(\frac{2}{r}-1\right) | \psi_0\rangle=12.\tag7
$$
Acting on just the right ket, as per your latter intention, also yields the same, 
$$
\langle \psi_0|\left( \frac{\delta(r)}{r^2} -\frac{2}{r} +\frac{2}{r^2} \right ) | \psi_0\rangle=12.\tag{7'}
$$
You never had to cancel the δ : it does not fail to know its place. It had better be, since the expectation of the commutator vanishes. 
But this helpful vanishing is not actually necessary for consistency. More generally, as per your off-diagonal example, 
$$
\langle \psi_1|\{p^2,1/r\}| \psi_0\rangle=  \langle \psi_1|\left(-\frac{1}{4}+ \frac{2}{r}\right)\frac{1}{r} + \frac{1}{r}\left(-1+\frac{2}{r}\right ) | \psi_0\rangle=86/27~. \tag8
$$
Alternatively, simply acting on the right, as done above, 
$$
\langle \psi_1|\left( \frac{\delta(r)}{r^2} -\frac{2}{r} +\frac{2}{r^2} \right) | \psi_0\rangle=\int dr e^{-3r/2}( \delta(r) -2r +2  )(2-r)  =86/27~.\tag{8'}
$$
Your two approaches are consistent, after all. The linchpin is that the  integration by parts implicit in the hermitean maneuver works free of surface terms at the origin.
This is not hard to generalize. For singular potentials (unlike yours, I understand), i.e. with $r^2V(r)$ not vanishing at the origin, take a look at 
Khelasvili & Nadareishvili 2010.
