First we need a vector space with an appropriate dot product. Consider the space of continuous functions on the interval [-1,1]. The dot product of f and g is the integral of f×g over the interval. If you like complex numbers, use f times the conjugate of g. We'll stick with the reals for now.
Let the powers of x act as a basis for the polynomial functions, then renormalize.
The first independent function is the constant function 1. There's nothing to change here.
The next independent function is x. The dot product of 1 and x is 0. That is, the integral of 1×x from -1 to 1 is 0. The functions are already orthogonal; there is nothing to do.
The next function is x2. This is already orthogonal to x. Subtract 1 × (1.x2/1.1), giving x2-(1/3).
The next function is x3, which is already orthogonal to the even degree Legendre polynomials. Subtract x × (x.x3/x.x) to get x3-(3/5)x.
The next function, derived from x4, becomes x4-(6/7)x2+(3/35). Continue this process, writing as many Legendre polynomials as you need.
Given a continuous function f on a closed interval, scale the domain so that the interval becomes [-1,1]. Take the dot product of f with each of the first n+1 Legendre polynomials. sometimes this integral can be computed directly; sometimes it is approximated by a computer. The result is the coefficient for that particular Legendre polynomial. Add these polynomials together to obtain a good approximation of f.
If g is the approximating polynomial of degree n, as computed above, let h = f-g, the error term. Let p be one of the first n+1 Legendre polynomials, and consider h.p. Remember that g is a linear combination of Legendre polynomials, all orthogonal. They all drop out, hence g.p is merely p.p times its coefficient. And that coefficient was computed from f.p, hence h.p = 0. The error term is orthogonal to the first n+1 Legendre polynomials, hence it is orthogonal to all polynomials of degree n or less.
Suppose we didn't really pick the best polynomial of degree n. Let g+e be an improvement. In other words, the integral of (h+e)2 is smaller than the integral of h2. Remember that an integral of a function squared is just that function dotted with itself. Expand (h+e).(h+e), remembering that h is orthogonal to every polynomial of degree n or less, hence h is orthogonal to e. The result is h.h+e.e. This is larger than h.h. Therefore g is the best nth degree polynomial approximation to f.