function [varargout]=orthzeroweights(n,a,b,c,win);

%  The function x=orthzeroweights(n,a,b,c,win) computes n roots of
%  the orthogonal poly. p_n(x) defined by thethree-term recurrence relation:
%      p_{k+1}(x)=(a(k+1) x -b(k+1)) p_k(x) -c(k+1) p_{k-1}(x), k>=0

%  The function [x,w]=orthzeroweights(n,a,b,c,win)  also returns the weights
%  of the associated Gauss-quadrature rule.
%
% Input: a,b,c: coefficients of the three term relation
%     win:  win=\int w(x) dx (where w(x) is the weight function)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:n
  dm(i)=b(i)/a(i);
end
for i=2:n
   ds(i-1)=sqrt(c(i)/(a(i-1)*a(i)));
end


J =diag(dm)+diag(ds,1)+diag(ds,-1);    %  Create Jacobi matrix
[e,x]= eig(J);               %  Compute eigenvalues
for i=1:n
dm(i)=x(i,i);
end
if nargout==1,
  varargout{1}=[dm];               %  Return n nodes
return;
end;
if  nargout==2,
  for i=1:n
    ds(i)=win*e(1,i)^2;
end
  varargout{1}=[dm];               %  Return n nodes
  varargout{2}=[ds];               %  Return n weights
  return;
end