restart;
# Find Gaussian Quadratue on [a,b] with weight function w
#
# Note: If the weight function is complicated, then the integrals 
# in the Gram-Schmidt algorithm may not have closed form solutions.
# In this case insert evalf( ... ) in suitable places to compute 
# the coeficients numerically.
with(LinearAlgebra):
n:=3;

a:=0;
b:=1;
w:=x->1;



for k from 0 to n
do
  q[k]:=x^k;
  for j from 0 to k-1
  do
    q[k]:=q[k]-q[j]*int(q[k]*q[j]*w(x),x=a..b);
  end:
  q[k]:=expand(q[k]/sqrt(int(q[k]^2*w(x),x=a..b)));
end:
p:=q[n];

c:=Vector([solve(p=0)]);

V:=VandermondeMatrix(c);

Y:=Vector([seq(int(x^k*w(x),x=a..b),k=0..n-1)]);

b:=LinearSolve(Transpose(V),Y);

F:=unapply(DotProduct(b,Map(f,Vector(c))),f);

F(x->1);
simplify(F(x->x));
simplify(F(x->x^2));
simplify(F(x->x^3));
simplify(F(x->x^4));
simplify(F(x->x^5));
simplify(F(x->x^6));