function xw=glquad(n)

% calculate zeros & weights for Gauss-Legendre quadrature
% by the Gollub-Welsch algorithm

A=zeros(n,n);
for i=2:n
   A(i-1,i)=(i-1)/sqrt((2*i-1)*(2*i-3));
end
A=A+A';

[V,D]=eig(A);

[x,in]=sort(diag(D));

for j=1:n,
    v=V(:,j);
    w(j)=2*v(1)^2;
end

w=w(in)';

xw=[x w];