function [v,xh,err] = poiss1d(N)
%  1D Poisson solve using DST - homogeneous Dirichlet BC
%
%  -u'' + u = f, u(0)=u(1)=0 
%
%  u(x) = exp(sin(pi*x))-1 
%  f(x) = pi^2(sin(pi*x)-cos(pi*x).^2).*exp(sin(pi*x))+exp(sin(pi*x))-1 

h = 1/N; 
xh = (0:N)'*h;
lam = 2*(1-cos(xh(2:N)*pi))/(h^2);

f = pi^2*(sin(pi*xh)-cos(pi*xh).^2).*exp(sin(pi*xh)) + ...
     (exp(sin(pi*xh))-1);

v = my_fftw_dst(f(2:N));  % transform   
v = v./(lam+1);     % Solve in Fourier space
v = my_fftw_dst(v);       % transform back
v = (1/(2*N))*[0; v; 0];      % set BC and normalize

x=0:0.01:1;  u=exp(sin(pi*x))-1;
err = norm(exp(sin(pi*xh))-1-v,'inf')
figure(1)
plot(x,u,xh,v,'o')
title([' N = ',num2str(N),'   error = ',num2str(err)]);


