function [nodes,weights]=gl_nodes_wts(a,b,n)

%     [nodes,weights]=gl_nodes_wts(a,b,n)
%
% produces the Gauss-Legendre nodes and weights 
% of order n on the interval [a,b].
%
% a,b     -  integration limits
% n       -  order of formula
% nodes   -  Gauss-Legendre nodes on [a,b]
% weights -  Gauss-Legendre weights on [a,b]
 
u=(1:n-1)./sqrt((2*(1:n-1)).^2-1);
[vc,bp]=eig(diag(u,-1)+diag(u,1));
[bp,k]=sort(diag(bp));
wf=2*vc(1,k)'.^2; 
nodes=(a+b)/2+((b-a)/2)*bp;
weights=(b-a)/2*wf;

end % gl_nodes_wts