The Legendre polynomial itself looks fine. To avoid the numerical imprecisions, you could work over QQ instead:

```
sage: legendre_P(90,t)(0.9)
1.56319477427853e13
sage: legendre_P(90,t)(QQ(0.9)).n()
-0.0874931768791928
```

It looks like you need 100+ bits of precision in order to evaluate the function in the original order:

```
sage: legendre_P(90, t)(RealIntervalField(100)("0.9"))
0.?e1
sage: legendre_P(90, t)(RealIntervalField(120)("0.9"))
-0.08749?
sage: legendre_P(90, t)(RealIntervalField(150)("0.9"))
-0.08749317687919?
```

The polynomial is very oscillatory in that there's a large amount of cancellation going on.

If you're really interested in the value, and only secondarily the polynomial, then you could use either `scipy`

or `mpmath`

, depending on how much precision you'll need:

```
sage: import scipy.special
sage: scipy.special.eval_legendre(90, 0.9)
-0.08749317687919235
sage: import mpmath
sage: mpmath.legendre(90, 0.9)
mpf('-0.08749317687919235')
```