How do we calculate the Lobatto gauss legendre points for a given order?

I have found an elegant code written to find the nodes of LGL polynomial extremely efficiently, but I am not able to understand the step where it uses Newton-Raphson method.

Can someone please help me understanding it? I tried to contact the author, but his email address is invalid.

 function [x,w,P]=lglnodes(N)

 % Truncation + 1
 N1=N+1;

 % Use the Chebyshev-Gauss-Lobatto nodes as the first guess
 x=cos(pi*(0:N)/N)';

 % The Legendre Vandermonde Matrix
 P=zeros(N1,N1);

 % Compute P_(N) using the recursion relation
 % Compute its first and second derivatives and 
 % update x using the Newton-Raphson method.

 xold=2;

 while max(abs(x-xold))>eps

     xold=x;

     P(:,1)=1;    P(:,2)=x;

     for k=2: N
         P(:,k+1)=( (2*k-1)*x.*P(:,k)-(k-1)*P(:,k-1) )/k;
     end

     x=xold-( x.*P(:,N1)-P(:,N) )./( N1*P(:,N1) );

 end

 w=2./(N*N1*P(:,N1).^2);

[Jason Turner]

I just posted about this on Math Stack Exchange, actually! I don't have an answer for you, I'm still working on that myself, but I was wondering if you had made any progress?