clear all
% -u''(x)+q(x)u(x)=f(x)
%  u(-1)=exp(-1);
%  u'(1)=exp(1)
q = @(x) x;
f = @(x) -exp(x)+x.*exp(x);
m = 10;
% Chebyshev + bordi
%x = [1,cos((2*(1:m-2)-1)/(2*(m-2))*pi),-1];
% Chebyshev-Lobatto
x = -cos((0:m-1)*pi/(m-1));
T(1,:) = ones(size(x))/sqrt(pi);
Tx(1,:) = zeros(size(x));
Txx(1,:) = zeros(size(x));
T(2,:) = x*sqrt(2/pi);
Tx(2,:) = ones(size(x))*sqrt(2/pi);
Txx(2,:) = zeros(size(x));
T(3,:) = 2*x.*T(2,:)-sqrt(2)*T(1,:);
Tx(3,:) = 2*T(2,:)+2*x.*Tx(2,:);
Txx(3,:) = 4*Tx(2,:);
for i = 4:m
  T(i,:) = 2*x.*T(i-1,:)-T(i-2,:);
  Tx(i,:) = 2*T(i-1,:)+2*x.*Tx(i-1,:)-Tx(i-2,:);
  Txx(i,:) = 2*Tx(i-1,:)+2*Tx(i-1,:)+2*x.*Txx(i-1,:)-Txx(i-2,:);
end
A = -Txx'+diag(q(x))*T';
b = f(x)';
A(1,:) = T(:,1)';
A(m,:) = Tx(:,m)';
b(1) = exp(-1);
b(m) = exp(1);
u = A\b;
norm(exp(x')-T'*u,inf)