It seems the special values of the Legendre Polynomials evaluated at zero are hard coded to be positive, but the sign should alternate. So, while legendre_P(0,0) correctly gives 1 and legendre_P(4,0) correctly gives 3/8, legendre_P(2,0) gives +1/2 instead of -1/2. Evaluating legendre_P(2,x) for any x close to zero gives a value close to -1/2 as it should. This problem is even evident on the documentation page (
https://paulmasson.github.io/sagemath-docs/functions/legendre_P.html) which shows the graph of P_2 passing through -1/2 but tabulates the special value at zero as +1/2.
Thanks.