function alfa = ChebyInterp(x,y)
% x and y are column n-vectors and x(1)<x(2)<...<x(n).
% alfa is a column n-vector with the property that if
% f(x) = alfa(1)B_{1}(x) +...+alfa(n)B_{n}(x)
% where the B_{1}(x),...,B_{n}(x) is the Chebyschev polynomial
% basis with respect to [a,b], then f(x(i)) = y(i), i=1:n.

n = length(x);
a = x(1);
b = x(n);
xtilde = -1+2*(x-a)/(b-a);
M = zeros(n,n);
M(:,1) = ones(n,1); % B_1(x)
M(:,2) = xtilde; % B_2(x)
for k=3:n
M(:,k) = 2*xtilde.*M(:,k-1)-M(:,k-2); % B_k(x)
end
% Solve the system M*alfa = y
[L,U,piv] = lu(M,'vector');
alfa = U\(L\y(piv));
% recursive but not vectorized -2
% column-row confusion -2
% lack of understanding of the problem (how the basis changes, or
% how the Cheb polys work) -3
% typo error but correct in general -1
% don't know how to use Lu -2