I am attempting to generate a as big as possible collection of orthogonal polynomials $p_1, p_2, ..., p_n$, $\left\langle p_i, q_i\right \rangle = \delta_{ij}$ where the inner product is with respect to a specific weight function.
It is well known that this can be done using the Gram-Schmidt process but my problem is that for large enough $n$ the resulting polynomials lose their orthogonality due to numerical instability. I have also implemented the Modified Gram-Schmidt algorithm but for around $n=30$ these also fail to be orthogonal.
My question is if there is a better method for generating these polynomials? I'm looking for a link or reference or description of the method.