restart;
dp:=(f,g)->int(f*g,x=-1..1);
nm:=f->sqrt(dp(f,f));
dp(x^2,x^4);
nm(x);
w[0]:=1;
v[0]:=w[0]/nm(w[0]);
w[1]:=x-dp(v[0],x)*v[0];
v[1]:=w[1]/nm(w[1]);
w[2]:=x^2-dp(v[0],x^2)*v[0]-dp(v[1],x^2)*v[1];
v[2]:=expand(w[2]/nm(w[2]),x);
N:=8;
for k from 0 to N do
   w[k]:=x^k;
   for j from 0 to k-1 do
       w[k]:=w[k]-dp(v[j],x^k)*v[j];
   end do;
   v[k]:=w[k]/nm(w[k]);
   print(k,v[k]);
end do:
plot({v[5]},x=-1..1);
R5:=solve(v[5]=0,x);
evalf(R5);
Digits:=16;
evalf(R5);