% NEWTCOTES Newton-Cotes formula
%
%    This generates the weights of a weighted n-point Newton-Cotes
%    quadrature formula with given nodes and weight function w. The
%    nodes are input through the column vector x with dimension >=n.
%    The two-column array ab of row-dimension >=floor((n+1)/2) is to 
%    hold the recursion coefficients of the orthogonal polynomials
%    with respect to the weight function w. The parameter eps0, close 
%    to machine precision, is used to detect a coincidence between the 
%    Gauss nodes and the Newton-Cotes nodes.
%
function w=NewtCotes(n,x,ab,eps0)
if size(ab,1)<floor((n+1)/2) | size(x)<n
  error('array x or ab too small')
  return
end
ng=floor((n+1)/2);
xw=gauss(ng,ab);
xg=xw(:,1); wg=xw(:,2);
if n==1
  w(1)=wg(1);
  return
end
lam=zeros(n,1); lam0=zeros(n,1); 
lam(1)=1;
for nu=2:n
  for mu=1:nu-1
    lam0(mu)=lam(mu);
  end
  for mu=1:nu-1
    lam(mu)=lam0(mu)/(x(mu)-x(nu));
  end
  lam(nu)=prod(1./(x(nu)-x(1:nu-1)));
end
w=zeros(n,1);
for nu=1:n
  for k=1:ng
    nu0=find(abs(xg(k)-x)<eps0);
    if isempty(nu0)
      al=lam(nu)/((xg(k)-x(nu))*sum(lam./(xg(k)-x)));
    elseif nu0==nu
      al=1;
    else
      al=0;
    end
    w(nu)=w(nu)+wg(k)*al;
  end
end
