function P3()
n = 16;
% Generate a problem based on n points
a = -2; b = 5; x = linspace(a,b,n)';
% Obtain and display the coefficients...
clc
format long
alfa = ChebyInterp(x,@MyF)
% Plot... 
m = 200; xvals = linspace(a,b,m)'; z = -1 + 2*(xvals-a)/(b-a);
Tkm2 = ones(m,1);
Tkm1 = z;
p = alfa(1)*Tkm2 + alfa(2)*Tkm1;
for k=3:n
    Tk = 2*z.*Tkm1 - Tkm2; p = p+alfa(k)*Tk; Tkm2 = Tkm1; Tkm1 = Tk;
end
plot(xvals,p,x,MyF(x),'o',xvals,MyF(xvals))
title(sprintf('n = %1d,  f(x) =  exp(-x/2)*(4*sin(2*x) - 3*sin(3*x))',n))
shg
 

function alfa = ChebyInterp(x,f);
% x is a column n-vector with a=x(1) < ... x(n)=b
% f is a handle to a function that is defined on the
% p(z) = alfa(1)T0(z) + ... +alfa(n) Tn-1(z) interpolates
% (x(1),f(x(1)),...,(x(n),f(x(n)))
n = length(x); a = x(1); b = x(n);
% Set up the linear systems...
A = zeros(n,n);
z = -1 + 2*(x-a)/(b-a);    % - 2 if you evaluated the T's at x instead of z
A(:,1) = ones(n,1);
A(:,2) = z;
for k=3:n
    A(:,k) = 2*z.*A(:,k-1) - A(:,k-2);
end
b = f(x);
% - 1 if your set-up involved O(n^3) flops. The above is O(n^2)
% Solve the linear system using LU...
[L,U,P] = lu(A);
alfa = U\(L\(P*b));
ConditionOfA = cond(A)

% test Function...
function y = MyF(x)
y = exp(-x/2).*(4*sin(2*x) - 3*sin(3*x));



