I'm trying to find the roots of Jacobi Polynomials using Roots[[JacobiP[n, α, β, x] == 0, x] (the roots are the abscissas for an n-point Jacobi-Gauss Quad rule). With n = 30, α = 10 and β = -0.95, I find that most of the roots are complex numbers and not real as I expected (correct me if I'm wrong to expect/assume that the roots of an n-order orthogonal polynomial such as a Jacobi polynomial are all real).
Using a MATLAB script that I found on the web, all the roots are real numbers for the same parameter values (n, α & β). Can anyone with experience in these issues provide an insight into the problem and a way to solve it if it is indeed a problem?
I'd really appreciate your input. Thank you in advance.