function spectralPa;
%
% Solution of the Poisson equation u'' + u = f(x) on (0,1) subject
% to homogeneous Neumann boundary conditions using second-order
% finite differences and spectral differentiation.
% The right-hand side is given by
% f(x) = (1.0+pi^2*sin(pi*x).^2-pi^2*cos(pi*x)) .* exp(cos(pi*x))
%

fprintf('\n');

for k = 1:10
    N = 2^k;
    [x,uFD] = SolveFD(N);
    [x,uSM] = SolveSM(N);
    uexact  = exp(cos(pi*x));
    Nk(k)   = N;
    aeFD(k) = max(abs(uFD-uexact));
    aeSM(k) = max(abs(uSM-uexact));
    fprintf(' k = %2d   error = %16.12e   error = %16.12e\n', ...
        k,aeFD(k),aeSM(k));
end;

figure;
plot(x,uSM);

figure;
loglog(Nk,aeFD,'b',Nk,aeSM,'r');

%
% Actual finite difference solver
%
function [x,u] = SolveFD(N);

h  = 1.0 / N;
x  = linspace(0.5*h,1-0.5*h,N)';
f = (1.0+pi^2*sin(pi*x).^2-pi^2*cos(pi*x)) .* exp(cos(pi*x));
DL = diag([-1,-2*ones(1,N-2),-1],0) + ...
     diag(ones(1,N-1),1) + diag(ones(1,N-1),-1);

A = DL / (h^2) + eye(N,N);
u = A  \ f;

%
% Actual spectral method solver
%
function [x,u] = SolveSM(N);

h    = 1.0 / N;
x    = linspace(0.5*h,1-0.5*h,N)';
f    = (1.0+pi^2*sin(pi*x).^2-pi^2*cos(pi*x)) .* exp(cos(pi*x));
ev   = -pi^2 * ((0:N-1)').^2 + 1.0;
fhat = dct(f);
uhat = fhat ./ ev;
u    = idct(uhat);