function [x,w]=gnodewt(n,a,b)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gaussian Integration points and coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[mu,v]=glnodewt(n);
w=0.5*(b-a)*v;
x=a+0.5*(mu+1)*(b-a);

function [x,w]=glnodewt(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gaussian Integration points and coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=(1:n-1)./sqrt(4*(1:n-1).^2-1);
J=diag(b,-1)+diag(b,1);
[V,D]=eig(J);[x,ix]=sort(diag(D));w=2*V(1,ix)'.^2;