function [x,w,varargout] = laguerrequad(n)
%
% Generate nodes and weights for
% Laguerre-Gauss quadrature.
% Note that x is a column vector
% and w is a row vector.
%
a = 2*(0:n-1)+1; % diagonal of J
u = 1:n-1; % upper diagonal of J
[V,Lambda] = eig(diag(u,1)+diag(a)+diag(u,-1));
[xx,i] = sort(diag(Lambda));
Vtop = V(1,:);
Vtop = Vtop(i);
ww = Vtop.^2;
%
i = find(ww>0);
w = ww(i);
x = xx(i);
nw = length(i);
if nw < n
   fprintf('\n *** Of %g Laguerre-Gauss weights, only %g are positive.\n\n',...
      n,nw)
end
if nargout==3
    varargout{1} = nw;
end