function lglnodes(N,a,b)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %
# % lglnodes.m
# %
# % Computes the Legendre-Gauss-Lobatto nodes, weights and the LGL Vandermonde 
# % matrix. The LGL nodes are the zeros of (1-x^2)*P'_N(x). Useful for numerical
# % integration and spectral methods. 
# %
# % Reference on LGL nodes and weights: 
# %   C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, "Spectral Methods
# %   in Fluid Dynamics," Section 2.3. Springer-Verlag 1987
# %
# % Written by Greg von Winckel - 04/17/2004
# % Contact: gregvw@chtm.unm.edu
# %
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#
# % Truncation + 1
	N1=N+1;

	#% Use the Chebyshev-Gauss-Lobatto nodes as the first guess
	x=cos.(pi*collect(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 maximum(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)/2 .* (b-a);
	x = (x .+ 1)*(b-a)/2.0 .+ a 
	return x,w,P
end