function x = LGL(N)

% Compute Legendre Gauss Lobatto points of order N
% Return in column vector x

Nl = N-1;
JN = zeros(Nl+1);

for j=2:Nl+1

  beta = (j-1)*(j+1)/((2*j+1)*(2*j-1));
  beta = sqrt(beta);
  
  JN(j,j-1) = beta;
  JN(j-1,j) = beta;
  
end;

[V,D] = eig(JN);

for j=1:Nl+1
  xp(j) = D(j,j);
end;

xh= sort(xp);

x = [-1.0,xh,1.0]';

return