function spectralPb;
%
% Solution of the Poisson equation u'' + u = 2*exp(x) on (0,1) subject
% to homogeneous Neumann boundary conditions using second-order
% finite differences and spectral differentiation
%

fprintf('\n');

for k = 1:10
    N = 2^k;
    [x,uFD] = SolveFD(N);
    [x,uSM] = SolveSM(N);
    uexact  = (exp(1)-cos(1))/sin(1) * cos(x) - sin(x) + exp(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  = 2.0*exp(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    = 2.0*exp(x);
ev   = -pi^2 * ((0:N-1)').^2 + 1.0;
fhat = dct(f);
uhat = fhat ./ ev;
u    = idct(uhat);