%% Example 6.12
% Calculation of Derivatives Using Chebyshev Derivative Matrix Operator


%% Initialization
close all;
clear N u j X U D j k

N  = 5;


u  = @(x) 4*x.^2.*(1-x.^2).*exp(-x/2);

j  = [0:N]';
X  = cos(j*pi/N);
U  = u(X);

D  = zeros(N+1,N+1);


%% Creation of the matrix operator

for j = [1:N+1]
    for k = [1:N+1]
        if j==k
            switch j
                case 1
                    D(j,k) = (2*N^2+1)/6;
                case N+1
                    D(j,k) = -(2*N^2+1)/6;
                otherwise
                    D(j,k) = -X(j)/(2*(1-X(j)^2));
            end;
        else
            if or(j==1,j==N+1) cj = 2; else cj = 1; end;
            if or(k==1,k==N+1) ck = 2; else ck = 1; end;
            D(j,k) = cj*(-1)^(j+k)/(ck*(X(j)-X(k)));
        end;
    end;
end;

%% Display the results

X
U
D
D*U