function [L,Lroots,GQw] = LegendrePolynomials(n)
%(Gram-Schmidt) L{k+1} = (x-b)*L{k} - c*L{k-1}
%(Gram-Schmidt) b=0 for Legendre Polynomials
%(Gram-Schmidt) L{k+1} = x*L{k} - c*L{k-1}
%(Gram-Schmidt) c = <x*L{k},L{k-1}> / <L{k-1},L{k-1}>>

L0 = 1;
L{1} = [1 0];
L{2} = [1 0 -1/3];

for k=1:2
  Lroots{k} = sort(roots(L{k}));
  %(scaling) We need <L{k},L{k}> = 2(2k+1) for the GQ-weight
  %(scaling) formula to be correct
  s = poly_int2(poly_mul(L{k},L{k})) * (2*k+1) / 2;
  for r=1:k
    xi = Lroots{k}(r);
    [Pprime,Pprimeval] = poly_diff(L{k},xi);
    GQw{k}(r) = 2 / (1-xi^2) / (Pprimeval^2) * s;
  end
  GQw{k} = reshape(GQw{k},length(GQw{k}),1);
end

for k = 2:(n-1)
  num = poly_mul(poly_mul([1 0],L{k}),L{k-1});
  den = poly_mul(L{k-1},L{k-1});
  c   = poly_int2(num) / poly_int2(den);
  L{k+1} = fliplr([L{k} 0]);
  nL  = length(L{k-1});
  L{k+1}(1:nL) = L{k+1}(1:nL) - c * fliplr(L{k-1});
  L{k+1} = fliplr(L{k+1});
  Lroots{k+1} = sort(roots(L{k+1}));
  s = poly_int2(poly_mul(L{k+1},L{k+1})) * (2*k+3) / 2;
  for r=1:(k+1)
    [Pprime,Pprimeval] = poly_diff(L{k+1},Lroots{k+1}(r));
    GQw{k+1}(r) = 2 / (1-Lroots{k+1}(r)^2) / (Pprimeval^2) * s;
  end
  GQw{k+1} = reshape(GQw{k+1},length(GQw{k+1}),1);
end 

for k=1:n
  L{k}   = reshape(L{k},length(L{k}),1);
end

return

function r = poly_mul(p,q)
pn = length(p);        qn = length(q);
p  = reshape(p,1,pn);  q  = reshape(q,1,qn);
p  = fliplr(p);        q  = fliplr(q);
r  = zeros(1,length(p)+length(q)-1);
for qk=1:qn
  r((1:pn)+qk-1) = r((1:pn)+qk-1) + p * q(qk);
end
r  = fliplr(r);
return

function q = poly_int(p)
pn = length(p);  p  = reshape(p,1,pn);  p  = fliplr(p);
q  = zeros(1,pn+1);
for k=1:pn
  q(k+1) = p(k)/(k);
end
q = fliplr(q);
return

function I = poly_int2(p)
q = poly_int(p);
I = polyval(q,1) - polyval(q,-1);
return

function [q,qx] = poly_diff(p,x)
pn = length(p);  p  = reshape(p,1,pn);  p  = fliplr(p);
q  = zeros(1,pn-1);
for k=1:(pn-1)
  q(k) = p(k+1)*k;
end
q = fliplr(q);
if (nargin>1 && nargout>1)
  qx = polyval(q,x);
end
return
