OK, this looks much better. Basically, the 'wrong' results was just floating point error, while the integral has really particular behavior near a=2. I'm not really a mathematician, but can anybody explain what goes on near a=2, i.e. how would you analyse the above integral as a function of 'a'.