@acer 

Thanks a lot for this. I've ended up using your first (method=_d01amc) approach.

The bug also continues through differentiation:

G:=(a)->int(exp(-(x^2-a)^2),x=0..infinity)

g:=(a)->D(G)(a)

evalf(g(0.54))
                         -0.3132963495

vs

F:=(a)->Int(exp(-(x^2-a)^2),x=0..infinity,method=_d01amc)

f:=(a)->D(F)(a)

evalf(f(0.54))
                        0.0009711319047

where f(a) is the correct value.

Oooh. That is a lot faster. What magic is the unapply doing?

@Preben Alsholm Thanks. That works for me.