%% Example 6.11
% Calculation of Integrals Using Discrete Chebyshev Transform


%% Initialization

close all;
clear N example u j X U a d I


format long;
N  = 16;
example = 'b';


% Functions we integrate
switch example
    case 'a'
        u  = @(y) (pi-1)/4 *  sin(0.5*(pi-1)*y + 0.5*(pi+1))./...
                              (0.5*(pi-1)*y + 0.5*(pi+1)).^3;

    case 'b'
        u  = @(y) 7/2 * log(3.5*y + 4.5)./(3.5*y + 4.5);
end;


%% Get the Chebyshev Transform
j  = [N:-1:0]';
X  = cos(j*pi/N);
U  = u(X);

a  = dct1(U);
a  = [a;0];
d  = zeros(N+2,1);

%% Numerical Integration
for j = [3:N+1]
    d(j) = (a(j-1)-a(j+1))/(2*(j-1));
end;
d(2)   = (2*a(1)-a(3))/2;
d(N+2) = a(N+1)/(N+1);

for j = [2:N+2]
    d(1) = d(1) + (-1)^j * d(j);
end;

%% I(x) is know directly in x = 1

I = sum(d)