function dc = demo_chebyshev()
    
    n    =  20;
    npts = 101;
    a    = -10;
    b    =  10;
    
    ci   = zeros(n);
    ti   = zeros(n);
    data = zeros(n,2);
    
    for i=1:n
        ci(i)=coeffc(i, n,a,b);
    end;
    
    step=(b-a)/(npts-1);

    x=a;
    for i=1:npts
        ti = arrayofTs( trans(x,a,b), n );
        y = 0.5*ci(1) + sum(ci(2:n).*ti(2:n));
        data(i,1) = x;
        data(i,2) = y;
        x=x+step;
    end;

    % fileID = fopen('vystup.dat','w');
    % fprintf(fileID,'%+18.14e    %+18.14e\n',data');
    % fclose(fileID);
    % data=load('vystup.dat');
   
    x1=(a:(b-a)/400:b);
    y1=func(x1);
    
    figure;
    plot(x1,y1,'b-');
    hold on;
    x2=data(:,1);
    y2=data(:,2);
    %figure(2);
    plot(x2,y2,'rx');
    hold off;
    data
    aa=axis;
    axis([a b aa(3) aa(4)]);
    
%%%%%%%%%%%%

function f=func(t)
% works for scalars as well as vectors  
    f = t.^4 .* log(10.0+abs(t)) .* cos(t);

    return;
    
%%%%%%%%%%%%

function tr=trans(x,a,b);

    tr = ( 2.0*x-(a+b) ) / (b-a);

    return;
    
%%%%%%%%%%%%

function tri=transinv(x,a,b)

    tri = (x*(b-a)+(a+b)) / 2.0;

    return;
    
%%%%%%%%%%%%

function c=coeffc(j, n,a,b);

    y=0.0;
    for i=1:n
        y = y + func( transinv( cos(pi/n*(i-0.5)), a, b) ) * cos(pi/n*(i-0.5)*(j-1));
    end
    c = 2.0*y/n;
    return;
    

%%%%%%%%%%%%

function Tvec = arrayofTs(x, n);

    Tvec(1)   = 1.0;
    Tvec(2)   = x;
    for i=3:n
        Tvec(i) = 2.0 * x * Tvec(i-1) - Tvec(i-2);
    end
    
    return;
    
%%%%%%%%%%%%
    
