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

h = 1/N;  xh = (0:N)'*h;
lam = -2*(1-cos(xh(1:N+1)*pi))/(h^2);
f =  pi^2*(-cos(pi*xh)+sin(pi*xh).^2).*exp(cos(pi*xh));
P = createcosP(N);

vh = P*f;                           % transform
vh(1) = 0;                          % set 0 mode  
vh(2:N+1) = vh(2:N+1)./lam(2:N+1);      
vh = P*vh;                          % transform back

vh = vh + (exp(1)-vh(1));
err = max(abs(vh - exp(cos(pi*xh))))




