function [x,w] = zplege(n)

% Gaussian-Legendre formula: x=Graussian points; w: integration weights
if( n <=1), disp( ' n must be > 1'); return;  end
[a,b] = coeflege(n);
JacM = diag(a) + diag(sqrt(b(2:n)),1) + diag(sqrt(b(2:n)),-1);
[w,x] = eig(JacM); x= diag(x); scal=2; w=w(1,:)'.^2*scal;
[x,ind]=sort(x); w=w(ind);