function [err,U]=bvpcheb(N);
% BVPCHEB  Solves linear BVP
%    u_xx + (sin x) u_x - 5 u = f(x),  u(-1)=1,  u_x(1)=0
% using a Chebyshev spectral method.  Here
%    f(x)= ((pi/2)^2+5) sin((pi/2) x) - (pi/2) sin(x) cos((pi/2) x)
% and the solution is known exactly:  u(x) = - sin((pi/2) x) 
% (it's a manufactured example).  
%   Requires CHEB to run.  This program is a modified version of P33 by
% Trefethen.  Note  "CHEB" and "P33" are from L. N. Trefethen, "Spectral 
% Methods in MATLAB", SIAM 2000.
%
% Format:
%   [err,U] = bvpcheb(N)
% where:
%     err = maximum error between exact and numerical soln
%       U = numerical solution at Cheb. points; need "x" from CHEB to plot
%       N = degree of Cheb. polynomial interpolation; N+1 points
%
% Example:
%    N=2:2:26; for j=1:13, err(j)=bvpcheb(N(j)); end
%    figure, semilogy(N,err,'o'),  grid on
%
% ELB 4/27/05
% See also:  BVPFD, CHEB, P33  

  [D,x] = cheb(N); D2 = D^2;     % D is (N+1) by (N+1); x is N+1 column
  D2(1,:) = D(1,:);              % Neumann condition at x = 1
  D2 = D2(1:N,1:N);              % D2 now N by N  
  D = D(1:N,1:N); D(1,:)=0;      % D modified to accomodate b.c.
  II = diag([0 ones(1,N-1)]);    % identity modified to accomodate b.c.
  % now produce matrix approx to operator  Lu = u_xx + (sin x) u_x - 5 u
  L = D2 + diag(sin(x(1:N))) * D - 5 * II;
  
  f=inline('((pi/2)^2+5)*sin((pi/2)*x)-(pi/2)*sin(x).*cos((pi/2)*x)','x');
  fp5 = 5 + f(x(2:N));      
  % v solves homogeneous prob: v_xx+(sin x)v_x-5v=5+f(x), v(-1)=0, v_x(1)=0                    
  v = L\[0; fp5];    % solve b.v.p. for v;  "0" here is "v_x(1)=0"
  U = 1+[v; 0];      % shift up to get soln u;  "0" here is "v(-1)=0"
  
  plot(x,U,'.','markersize',18)
  xx = -1:.01:1;  
  % next method of computing smooth version of the numerical solution is
  % not so good if N is larger than 25 or so;  for plotting use alternate
  uu = polyval(polyfit(x,U,N),xx);    
  % uu = interp1(x,u,xx,'cubic');   % alternate for *plotting* if N > 25
  hold on,  plot(xx,uu),  grid on
  uexact = -sin((pi/2)*xx); 
  plot(xx,uexact,'r'),  hold off
  
  err=norm(uu-uexact,inf);  
  % err=norm(U+sin((pi/2)*x),'inf');  % alt. error computation if N > 25
  title(['max err = ' num2str(err)],'fontsize',12)