(* Gram-Schmidt Orthogonalization to get the orthogonal polynomials *)
a = -1; b = 1;

w[x_] := Sqrt[1-x^2]

norm[n_] := π/2

(* These are the inputs: need to be changed to get the different polynomials.
a,b integration limit, w weight function and norm gives the normalization *)

LP[0, x_] := Sqrt[norm[0]]/Sqrt[Integrate[w[x], {x, a, b}]]

For[n=1, n ≤ 5, n++, 
GS[n_, x_] := Sum[Integrate[x^n LP[i, x] w[x], {x, a, b}] LP[i, x] / norm[i], {i, 0, n-1}];
LPP[n_, x_] := x^n - GS[n, x];
LP[n_, x_] := LPP[n, x] Sqrt[norm[n]] / Sqrt[Integrate[LPP[n, x]^2 w[x], {x, a, b}]]]

MatrixForm[Table[{Simplify[LP[n, x]], ChebyshevU[n, x]}, {n, 0, 5}]]

(* Fun with Fourier expansion *)

(* We try Fourier expansion for f(x) = 0, -1 ≤ x < 0, f(x) = 1, 0 ≤ x ≤ 1. This is a sine series for the odd function f(x) - 0.5. *)

f[x_] := 0.5

For[n= 1, n ≤ 10, n++,
fs[n_] := FourierSinCoefficient[f[x], x, n, FourierParameters -> {1, Pi}] ]

Table[{fs[n]}, {n, 1, 10}]

fj[x_] := 0.5 + Sum[fs[n] Sin[n Pi x], {n, 1, 10}]

Plot[Evaluate[fj[x]], {x, 0, 1}]

fl[x_, m_] := 4 Sum[Sin[n Pi x] / n, {n, 1, m, 2}]

Plot[{fl[x, 9]/Pi, fl[x, 51]/Pi}, {x, 0, 1}, PlotRange -> {0,2}]

fp[x_, n_] := 0.5+FourierSinSeries[0.5, x, n, FourierParameters -> {1, Pi}]

Plot[Evaluate[fp[x, 5]], {x, -2, 2}]

Plot[Evaluate[FourierSinSeries[1, x, 5, FourierParameters -> {1, Pi}]], {x, 0, 1}]

Table[NIntegrate[(1-fp[x,n])^2, {x, 0, 1}], {n, 10, 50, 20}]

(* If we wanted to write the full function as defined in the problem, the notation is as follows. *)

g[x_] := Piecewise[{{0, -1<x<0},{1, 0 ≤ x ≤ 1}}]

Plot[g[x], {x, -1, 1}]

(* Could also expand the function in Legendre basis. *)

nmax=20

lp[n_] := NIntegrate[LegendreP[n, x], {x, 0, 1}]

Table[lp[n], {n, 1, nmax, 2}]

lg1[x_] := 0.5+Sum[0.5*(2 n+1)*lp[n]*LegendreP[n, x], {n, 1, nmax, 2}]

Plot[Evaluate[lg1[x]], {x, -1, 1}]