function [x,u]=ChebyshevPoisson1d(f,uL,uR,n)
% CHEBYSHEVPOISSON1D Chebyshev spectral solver for the 1d Poisson equation 
%   [x,u]=ChebyshevPoisson1d(f,uL,uR,n) solves uxx=f on [-1,1] with
%   u(-1)=uL, u(1)=uR using n points (excluding end points)

n=n+1;                                      % Number of intervals
xi=(0:pi/n:pi)';                            % Interpolation points
ff=f(cos(xi));                              % Sample f at Chebyshev pts
fhat=fft([ff(1:end-1);ff(end:-1:2)])/(2*n); % Discrete Cosine Transform
fhat=fhat(1:n-1);                           % (first n-1 coeffs only)
fhat(2:end)=fhat(2:end)*2;

D=zeros(n+1);                               % Differentiation matrix
for j=1:n,
  D(j:-2:1,j+1)=2*j;
end;
D(1,:)=D(1,:)/2;
A=D^2;                                      % Square to get Laplacian
A(end-1,:)=(-1).^(0:n);                     % Left boundary condition
A(end,:)=ones(1,n+1);                       % Right boundary condition
b=[fhat; uL; uR];                           % Right hand side
uhat=A\b;
x=cos((0:pi/n:pi)');                        % Sample points
uuhat=0.5*[uhat(1:end-1);uhat(end:-1:2)];   % Inverse DCT
uuhat(1)=uuhat(1)*2; uuhat(n+1)=uuhat(n+1)*2;
u=ifft(uuhat)*2*n;
u=u(1:n+1);