% Solve linear BVP 
%     q'' = cos(pi*x/2), -1 < x < 1
%     q'(-1)-q(-1)=0, q'(1)+q(1)=0
% with Chebyshev PS method of orders N = 2^p, p = 1,...,4

  prange = 1:4;
  xx = -1:0.01:1;
  qq = (2/pi)*((2/pi)*cos(pi*xx/2)+1);
  err = zeros(1,length(prange));
  errfd = zeros(1,length(prange));

  for p = prange
    N = 2^p;
    [D,x] = cheb(N);
    D2 = D^2;              
    D2(1,:) = D(1,:);   
    D2(1,1) = D2(1,1)+1;              % D2(1,:)q = q'+q at x=1
    D2(N+1,:) = D(N+1,:);    
    D2(N+1,N+1) = D2(N+1,N+1)-1;      % D2(N+1,:)q = q'-q at x=-1

    D2fd = D2; 
    
    % For FD solution, replace spectral BCs with FD BCs.
    D2fd(1,:) = 0;
    D2fd(1,1) = 1/(x(1) - x(2)) + 1/(x(1)-x(3));
    D2fd(1,1) = D2fd(1,1)+1;
    D2fd(1,2) = (x(1) - x(3))/((x(2) - x(1))*(x(2) - x(3)));
    D2fd(1,3) = (x(1) - x(2))/((x(3) - x(1))*(x(3) - x(2)));
    D2fd(N+1,:) = 0;
    D2fd(N+1,N+1) = 1/(x(N+1) - x(N)) + 1/(x(N+1)-x(N-1));
    D2fd(N+1,N+1) = D2fd(N+1,N+1)-1;
    D2fd(N+1,N) = (x(N+1) - x(N-1))/((x(N) - x(N+1))*(x(N) - x(N-1)));
    D2fd(N+1,N-1) = (x(N+1) - x(N))/((x(N-1) - x(N+1))*(x(N-1) - x(N)));

    f = [0; -cos(pi*x(2:N)/2); 0];      % RHS -cos(pi*x/2) except at endpoints,
                                      % where f is dictated by BCs.
    q = D2\f;                         % Numerical soln with Chebyshev BCs.
    qfd = D2fd\f;                     % Numerical soln with FD BCs
    qex = (2/pi)*((2/pi)*cos(pi*x/2)+1);  % Exact soln.
    err(p) = norm(q-qex,inf);
    errfd(p) = norm(qfd-qex,inf);
  end

  % Plot of Chebyshev solution with Chebyshev BCs
  
  subplot(2,2,1)
  xlabel('x')
  ylabel('q')
  plot(x,q,'r.',xx,qq,'b-')
  legend('Chebyshev','Exact')
  title(['N=' int2str(N)])
  
  % Comparison of max-norm error convergence with Chebyshev and FD
  % BCs.
  
  subplot(2,2,2)
  np = length(prange);
  Nrange = 2.^prange;
  beta = -(log(err(np))-log(err(np-1)))/(Nrange(np) - Nrange(np-1));
  errfit = err(np)*exp(beta*(Nrange(np) - Nrange));
  errfdfit = errfd(np)*(Nrange(np)./Nrange).^4;
  semilogy(Nrange,err,'rx',Nrange,errfit,'r--',Nrange,errfd,'b+',Nrange,errfdfit,'b-')
  xlabel('N')
  ylabel('Max norm error')
  ylim([1e-20, 10])
  legend('Spectral BCs','b e^{-2.15N} fit','FD BCs','cN^{-4} fit',3)