@mmcdara Thanks for your reply, it contains a lot of interesting knowledge I can use to improve the resolution of the problem. However, I have to note some things:
First point: I don't know why you assume anything about my definition of HN, indeed, it is perfectly well defined. HN are the normalized Hermite functions, which are built from the Hermite polynomials and their inner product. You can easily check that there is no such "unjustified sqrt": https://en.wikipedia.org/wiki/Hermite_polynomials#Hermite_functions
Second point: If you read carefully the answer you refer to, you will find the sentence "Regarding the RootOf function, fsolve would work nice in that case. I will change that". I agree with the usage of fsolve for that task, and I have already added that change to my to-do list.
Therefore, my "pure nonsense" answer about complex solutions is obviously not about finding the roots of the polynomials for any "complex gaussian quadrature"... It refers to the solutions of the final nonlinear system for the coefficients of the pseudospectral expansion, which can be, and indeed some of they are, complex. Those solutions containing non-zero imaginary parts are then discarded from the process of computing the pseudospectral expansion.