function [x,w,varargout] = legendrequad01(n,varargin)
%
% If [x,w] = legendrequad01(n), and f is a
% function defined on [0,1], then the formula
% w*f(x) approximates the integral of f over
% the interval [0,1].  Here w is a row vector
% of weights and x is a column vector of nodes
% for Legendre-Gauss quadrature on [0,1].
%
% If the function f is defined on an interval
% [alpha,beta], you can use the expression
% (beta-alpha)*w*f((beta-alpha)*x+alpha)
% to approximate the integral of f over
% [alpha,beta].  This amounts to doing a change
% of variable to convert the integral over
% [alpha,beta] into an integral over [0,1].
% Explicitly coding this formula is most
% efficient when you have to integrate a
% function (or functions) over several
% different intervals, because you only have
% to generate the nodes and weights one time.
%
% However, if you just need to do one integral
% over an interval [alpha,beta], you can use
% [x,w] = legendrequad01(n,alpha,beta)
% and then w*f(x) to approximate the integral
% of f over [alpha,beta].  Providing alpha
% and beta to legendrequad01 tells it to
% replace w with (beta-alpha)*w and to replace
% x with (beta-alpha)*x+alpha.

a = repmat(1/2,1,n);                % main diagonal of J
u = sqrt(1./(4*(4-1./(1:n-1).^2))); % 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 Legendre-Gauss weights on [0,1], only %g are positive.\n\n',...
      n,nw)
end
if nargout==3
    varargout{1} = nw;
end

if nargin>1               % added July 30, 2011.
   alpha = varargin{1};
   beta = varargin{2};
   bma = beta-alpha;
   w = bma*w;
   x = bma*x+alpha;
end