I solved it on paper using the trick mentioned above. Given that h() is a product of three L2 functions, we can use Parseval's on the first two + last, and the Fourier convolution on the first two. The new integrand showed that it is zero everywhere except for m=p and q=n or q=n+1.

FInally, the integral is I/2/Pi/m for m=p, q=n; and -I/2/Pi/m for m=p, q=n+1 ... and 0 otherwise. This matches numerical integration.

I attached the MW file with my derivations and some explanations for anyone curious (most expressions were derived on paper and input in Maple by hand)

psi-phi-phi-2.mw

I still would like to hear the opinion of the Maple experts as to why Maple gives inconsistent results for the initial h(x) integral ...

It also gets stuck for the equivalent integral of Fourier transforms, which is why I had to resort to pen and paper for pretty much all steps.

Alex.