function xp = ChebyshevNodes(xstart, xend, npoints)

    xp = zeros(npoints, 1) ;
    length = xend-xstart ;

    for ipoint = 1 : npoints
        xp(ipoint) = -length*(cos(pi*(ipoint-1)/(npoints-1))+1)/2+xstart ;
    end
      
end

function [PHI, PHIx, PHIxx] = ChebyshevInterpolation(npoints, xquery)
  
    nquery = numel(xquery) ;
    length = max(xquery)-min(xquery) ;
    
    xq = ChebyshevNodes(-length/2, +length/2, nquery) ;

    PHI = zeros(nquery, npoints) ;
    PHIx = zeros(nquery, npoints) ;
    PHIxx = zeros(nquery, npoints) ;

    PHI(:, 1) = 1 ; PHI(:, 2) = 2*xq/length ;
    PHIx(:, 1) = 0 ; PHIx(:, 2) = 2/length ;
    PHIxx(:, 1) = 0 ; PHIxx(:, 2) = 0 ;

    for ipoint = 3 : npoints
        PHI(:, ipoint) = 2*PHI(:, 2).*PHI(:, ipoint-1)-PHI(:, ipoint-2) ;
        PHIx(:, ipoint) = 4*PHI(:, ipoint-1)/length ...
            + 2*PHI(:, 2).*PHIx(:, ipoint-1)-PHIx(:, ipoint-2) ;
        PHIxx(:, ipoint) = 8*PHIx(:, ipoint-1)/length ...
            + 2*PHI(:, 2).*PHIxx(:, ipoint-1)-PHIxx(:, ipoint-2) ;
    end
    
end