function c=spectbvp(K,dj,xend);
% Uses orthogonal fcns based on decaying exponentials and collocation
% to solve boundary value problem stated in A#6 Problem IX.
% Calling:     c=spectbvp(K,dj,xend)     uses K exponentials 
% exp(-dj x),...,exp(-K dj x) and computes on interval [0 xend]
% for finite difference comparison.  Returns coeffs of soln.
% Method described in Problem IX:   >> spectbvp(5,1,3)
% Also try:  >> spectbvp(20,.25,5)     ELB 11/8/03

% Compute R by method due to Dixon Jones:
[k j]=meshgrid(1:K,1:K); R0=2*gammaln(k)-gammaln(abs(k-j)+1)-gammaln(k+j+1);
RK=triu(exp(R0+.5*log(2*k.*k.*j)))/sqrt(dj);
% Alternative fails to get R when K>12 or so because chol fails:
% B=hilb(K+1)/dj; B=B(1:K,2:K+1); RK=chol(B);  

% roots z of qK(x) are collocation points:
emz=roots(flipud(RK\([zeros(1,K-1) 1]'))); z=-log(emz)/dj;
R=RK(1:K-1,1:K-1);  % only need subblock of R from now on
Aco=zeros(K-1,K-1); for j=1:K-1, Aco(:,j)=exp(-j*dj*z); end
Qco=Aco/R; % q_j at colloc pts
ddQ=(Aco*diag([(dj*(1:K-1)).^2]))/R; % q_j'' at colloc pts
M=[ones(1,K-1)/R;  -ddQ+repmat(5+4*cos(z),1,K-1).*Qco];  
% M is overdetermined sys (1) and (2); remove last equation, solve
c=M(1:K-1,:)\([1 zeros(1,K-2)]');  % c are coeffs of q_j to form uhat

% for plot of soln, discretize x variable
n=801; x=(linspace(0,xend,n))';
A=zeros(n,K-1); for j=1:K-1,  A(:,j)=exp(-j*dj*x); end % the exponentials
Q=A/R;  uhat=Q*c;  % uhat is spectral solution as a fcn of x

% now construct finite difference solns and compare
for N=[5 23]
    dxx=xend/N; xx=(dxx:dxx:xend-dxx)'; % xx is loc of interior pts
    e = ones(N-1,1); DD = spdiags([e -2*e e], -1:1, N-1,N-1);
    DD=-(1/(dxx*dxx))*DD + spdiags(5+4*cos(xx),0,N-1,N-1);
    if N==5, xxa=xx; UUa=DD\([1/(dxx*dxx) zeros(1,N-2)]');
    else, xxb=xx; UUb=DD\([1/(dxx*dxx) zeros(1,N-2)]'); end
end
subplot(2,1,1), semilogy(x,uhat,xxa,UUa,'go',xxb,UUb,'ro'), hold on
legend('spectral soln','finite difference soln','finer finite diff. soln')
semilogy(x,exp(-x),'--',x,exp(-3*x),'--')
text(xend/2,1.5*exp(-xend/2),'e^{-x}')
text(xend/2,.5*exp(-3*xend/2),'e^{-3x}')
title(['Three solns to a boundary value problem.  '...
    'Seen in log scale, spectral soln has no "end effect"'])
hold off, subplot(2,1,2) 
diffu=-(A*diag([(dj*(1:K-1)).^2]))/R+repmat(5+4*cos(x),1,K-1).*Q;
plot(x,diffu*c), grid on
title('Residual   -u" + (5 + 4 cos x) u   for spectral soln.')