I guess this question has sat here awhile. I just stumbled on it during a search. The reason is that adaptive quadrature works by sampling function values at a finite number of points. The fundamental assumption is that the initial sampling will reveal at least some of the function variation. If the initial points all return a function value of zero, or near to it, then the integrator will assume that the function is a constant zero function (or close enough) and return an integral of zero.
There are a couple of ways out of this. The brute force and ignorance way is to divide your region up into subsets with similar amounts of variation, making sure that you the integrator can't help but notice the function variation in the areas where there is a lot going on. Then down in the "tail" of the function, where the function values are much smaller, returning zero is OK. The other way is to transform the infinite interval into a finite one. This results in at least a mild singularity, but some integrators can handle it. DBLQUAD doesn't really. The new INTEGRAL2 function can often handle improper integrals. Just try putting -inf and inf in there instead of large finite numbers.
