      function [z,w] = zwgll(N);

% Gauss-Lobatto-Legendre points z and weights w.

%    z are the N+1 zeros of  (1-z^2)*P'_N(z)
%    w are the weights.


n = N+1; w(1:n)=0; z(1:n)=0; z(1)=-1; z(n)= 1;

if N>1,
 if N==2 
   z(2)=0; 
 else
  M=zeros(N-1,N-1);
  for i=1:N-2,
    M(i,i+1)=(1/2)*sqrt((i*(i+2))/((i+1/2)*(i+3/2)));
    M(i+1,i)=M(i,i+1);
  end;
  z(2:N)=sort(eig(M));  % Eigenvalues of companion matrix
 end;
end;

%compute the weights w

w(1)=2/(N*n); w(n)=w(1);

for i=2:N, x=z(i);
  z0=1; z1=x;
  for j=1:N-1,
    z2=x.*z1*(2*j+1)/(j+1)-z0*j/(j+1);
    z0=z1; z1=z2;
  end;
  w(i)=2/(N*(n)*z2*z2);
end;

z=z'; w=w'; 

