@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.

@acer One solution is not enough, that's why I tried with solve. I need to get all the possible solutions, then discard complex ones and study the behavior of the ones that only provide pure real coefficients for the expansion.

Regarding the RootOf function, fsolve would work nice in that case. I will change that.

@Carl Love Yes, sorry. When running the solve function, I get a "lost kernel connection" message after several seconds of execution. I read that this might be a problem related to a corrupted install, but I reinstalled Maple and the problem persists.

@Mariusz Iwaniuk I get the same output using _CCquad. However, I realised I had

Digits := 50:

at the beginning of the whole worksheet. Apparently, removing this line solves the problem. I used that variable to increase the resolution of Maple, so how can I get such resolution and still get the proper output from the numerical integration?

@C_R I'm working with 2021.1, could that be the reason?