When computing the dot product of a Bra and a Ket from two different basis we are
implicitly referring to an object operating in the tensor product space. If we want to
compute

Everything will work if the product of the bra from one basis with the ket of an other
basis is defined as



then we would have




That is

After testing the submited intdeo integration function I conclude that this is the way to go.
Using Ooura BK0 approximation it is possible to speedup the calculation by an average
factor of 20 without lossing precision.

I would vote for Axel's solution but since my reputation in 0, I am not allowed to vote

Thank you again.

1) Dont use maxintervals, you will eat all your memory (onless you have only one value to compute).
    Split the integration interval
2) The results for this case appears good with a relative difference of 1.2e-10.

The kmin=0.000224399475256414 to kmax =91.5632886852008 was splitted in 52 intervals.
NAG (_d01akc) + Ooura's BK0 now gives -0.000109403613448212
     compared with the previous one        -0.000109403613457678  obtained with maxintervals=1000000

putting kl=kmax in the following expression for the tail of the integral form kmax to infinity
with the more precise values for 

rho1=0.323147799002253 and a = 0.01785 and dz=7000

I obtain NAG (_d01akc) + Ooura's BK0 kmin..infinity = −0.000109403613448236

Axel Vogt used used a=0.323 and obtained  intdeo + BesselK(0,x) : -0.000109408897877997
For this case Ooura BK0 + NAG(_d01akc) gives                                -0.000109408897865057

This appears coherent.

Again thanks a lot for your help.

I think I have the solution to my problem to this unusual memory consumption.
This is not elegant in the mathematical sense but practical.

in evalf(Int(f, kmin..kmax, method = _d01akc,epsilon=eps,maxintervals=1000000))

maxintervals=1000000 appears to create the problem. I will never use it again.

I used

evalf(Int(f, kmin..kmax, method = _d01akc,epsilon=eps) # eps =1e-12

and divided the interval kmin..kmax in subintervals not presenting more than 2000 oscillations
and avoiding this way the default maximum subinterval limit of 500 in _d01akc.



Thank you for your support and good advice.

Im is the imaginary part of complex number.
Im (A+B*I) = B assuming A and B real.