Orthogonal Polynomials: Mathematica

Orthogonal polynomials are everywhere.

The following Mathematica program takes in a weight function (as a function of x), the domain (a and b), and spits out the first "n" corresponding orthonormal polynomials.

OrthoPoly[w_, a_, b_, n_] := Module[{monoBasis, T},
  monoBasis = Table[{x^i}, {i, 0, n - 1}];
  oP = monoBasis;
  oP[[1]] = oP[[1]]/Sqrt[Integrate[w*oP[[1]]*oP[[1]], {x, a, b}]];
  For[i = 2, i <= n, i = i + 1,
   For[T = 0; j = 1, j < i, j = j + 1,
    T = T + Integrate[w*oP[[i]]*oP[[j]], {x, a, b}]*oP[[j]];
    ];
   oP[[i]] = oP[[i]] - T;
   oP[[i]] =
    oP[[i]]/Sqrt[Integrate[w*oP[[i]]*oP[[i]], {x, a, b}]] //
     Simplify;
   ] ;
  oP
  ]

Here's a screenshot, for the first four Chebyshev and Legendre polynomials (click to enlarge). Note that these polynomials are unique up to a multiplicative constant. Orthonormality freezes that constant.