   
function [D, x] = LegDifferenceMatrix(N)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % DifferenceMatrix.m
    %
    % Computes the Legendre differentiation matrix with collocation at the 
    % Legendre-Gauss-Lobatto nodes.
    %
    % Reference: 
    %   C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Tang, "Spectral Methods
    %   in Fluid Dynamics," Section 2.3. Springer-Verlag 1987
    %
    % Written by Houwang Tu - 01/22/2021
    % Contact: tuhouwang96@163.com
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    [x, ~, ~] = LGLnodes(N);
    
    % The i-th row and j-th column of L and LP stores L_j(x_i)    
    [L, ~]    = legpoly(N, x);
    D  = zeros(N + 1, N + 1);

    for k = 1 : N + 1
        for j = 1 : N + 1
            if(k ~= j)
                D(k, j) = L(k, N + 1) / L(j, N + 1) / (x(k) - x(j));
            end
        end
    end
    
    D(1, 1)         =   N * (N + 1) / 4;
    D(N + 1, N + 1) = - N * (N + 1) / 4;
end
 


%{
function [D, x]=DifferenceMatrix(N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% legDc.m
%
% Computes the Legendre differentiation matrix with collocation at the 
% Legendre-Gauss-Lobatto nodes.
%
% Reference: 
%   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 - 05/26/2004
% Contact: gregvw@chtm.unm.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Truncation + 1
N1=N+1;
% CGL nodes
xc=cos(pi*(0:N)/N)';
% Uniform nodes
xu=linspace(-1,1,N1)';
% Make a close first guess to reduce iterations
if N<3
    x=xc;
else
    x=xc+sin(pi*xu)./(4*N);
end
% 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
X=repmat(x,1,N1);
Xdiff=X-X'+eye(N1);
L=repmat(P(:,N1),1,N1);
L(1:(N1+1):N1*N1)=1;
D=(L./(Xdiff.*L'));
D(1:(N1+1):N1*N1)=0;
D(1)=(N1*N)/4;
D(N1*N1)=-(N1*N)/4;
    
end
%}
