For the Legendre polynomials orthogonality requires the weighting function x = 1. I think that you are confusing the Legendre polynomials for the Chebyshev polynomial of the first kind that are orthogonal with the weight function 1/sqrt(1-x^2).
Although the inner product test of orthogonality may seem crude, it is practical because it is how the polynomials will be used. The ideal result should converge as the number of points increase and as one can see this case does not occur.
I also coded my own legendre dissociated function whose plots match the wiki article that I cited above. It still does not achieve the orthogonal results that I am seeking.
The function is below.
function y = LegPoly(x, df)
[m, n] = size(x); if m < n x = x'; m = n; end
y = zeros(m, df+1); y(:, 1) = ones(m, 1); y(:, 2) = x;
for i = 3:df+1 j = i-1; y(:, i) = ((2*j+1).*x.*y(:, j) - j*y(:, j-1))/i; end