Clear[ips,P,a,b];
ips[f_,g_] := Integrate[f g, {x,-1,1}];
a[n_] := ips[ x P[n-1,x] , P[n-1,x] ] / ips[ P[n-1,x], P[n-1,x] ]; 
b[n_] := ips[ x P[n-1,x] , P[n-2,x] ] / ips[ P[n-2,x], P[n-2,x] ];
P[0,x_] = 1; 
P[1,x_] = x - a[1];
P[n_,x_] := P[n,x] = (x-a[n]) P[n-1,x] - b[n] P[n-2,x]   ;

For[j=0, j<=5, j++,
   Print["P",j,"(x)=" Simplify[P[j,x]]] ]

For[j=0, j<=5, j++,
   Print["P",j,"(x)=",  LegendreP[j,x]] ]

BestPoly[f_,p_] := Simplify[
        Sum[ (ips[ f[x] , P[k,x] ] /
        ips[ P[k,x], P[k,x] ]) P[k,x]  ,{k,0,p} ] ]

Clear[f,g];
f[x_] = Exp[x];
g[x_]=BestPoly[f,2] //N

Clear[func,approx];
func=Plot[f[x],{x,-1,1}, DisplayFunction->Identity]; 
approx=Plot[g[x],{x,-1,1}, DisplayFunction->Identity];
Show[func,approx,DisplayFunction->$DisplayFunction]

Clear[f,g];
f[x_] = Cos[4 x]; 
g[x_]=BestPoly[f,4] //N

Clear[func,approx];
func=Plot[f[x],{x,-1,1}, DisplayFunction->Identity];
approx=Plot[g[x],{x,-1,1}, DisplayFunction->Identity];
Show[func,approx,DisplayFunction->$DisplayFunction]