I've reported this at http://trac.sagemath.org/ticket/18059. I don't have an answer (yet) for why this is happening, because ordinarily it indicates there is something like a complex number as an intermediate result. But I will note that numerical_integral(lambda r: y.subs(theta=1,phi=1),0,5) yields this problem while z = y.subs(theta=pi/2,phi=pi/2); numerical_integral(z,0,5) does not, and yet numerical_integral(lambda r: z,0,5) does but numerical_integral(lambda r: 0.240000000000000*e^(-1/20*r)/(4*pi - pi^2 + 2*r^2 + 8), 0,5) is fine!