As cyclist has noted, the factor 1/(1-x^2) is essential in performing an orthogonality test. The Associated Legendre "polynomials" for differing m values are only orthogonal when each function is divided by sqrt(1-x^2).
Also note that when you make this change, your y*y' approximation to an integral would then yield NaNs at the two endpoints because of a zero-divided-by-zero occurrence. In the case of m1 = 0 and m2 = 1, you will get actual square root singularities in the integrand at the two ends even though the integral is still finite. To get accurate results for the orthogonality test I would advise you to use for-loops to compute with matlab's numerical quadrature function, 'integral', which can handle singularities, rather than your crude y*y' method. Either that or make an appropriate change of variable so that points are packed closely together at the two endpoints.