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]