function xw=gauss(N,ab)
%---------------------------------------------------------------------- 
% This function computes nodes and weights of the gaussian quadrature formula.
%
% INPUTS:
% N - number of nodes
% ab - n alpha coefficients and n beta coefficients of monic Jacobi polynomials
%
% OUTPUTS:
% xw - matrix with two columns: in the first column there are the nodes of gaussian
% quadrature formula, in the second one there are the quadrature weights
% 
% References:
% Walter Gautschi website


N0=size(ab,1); if N0<N, error('input array ab too short'), end
J=zeros(N);
for n=1:N, J(n,n)=ab(n,1); end
for n=2:N
  J(n,n-1)=sqrt(ab(n,2));
  J(n-1,n)=J(n,n-1);
end
[V,D]=eig(J);
[D,I]=sort(diag(D));
V=V(:,I);
xw=[D ab(1,2)*V(1,:)'.^2];
