High precision evaluation of orthogonal polynomials

Hi everybody,

I need to evaluate some orthogonal polynomials with high precision. For example for Legendre polynomials this is what I do:

 % psi=@(n,x)(feval(symengine,'orthpoly::legendre',n,x));

But the problem is that the function does not distinguish between inputs which are different in digit numbers beyond double precision. For example psi(1,0.333333333333333324234) and psi(1,0.33333333333333345) have the same answer, while in my calculations I need different answers which are more precise.

Any suggestions?!

Thanks very much in advance!

[Mehdi]

I think the problem is that at first, an input is converted to a fraction which is not different for inputs which are different beyond double precision digits...