function [vh,xh] = poiss1d_4P(f,N)
%
%  2nd order 1D Poisson solve with DST 
%  implemented by matrix multiply by P 
%
%  u'' = f,  u(0)=u(1)=0 
%

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

h = 1/N;          
xh = (0:N)'*h;               % discrete grid
fh = f(xh);
lam = -2*(1-cos(xh(2:N)*pi))/(h^2);  % evals 
P = createP(N); 

fh = fh(2:N)+(1/12)*(fh(1:N-1)-2*fh(2:N)+fh(3:N+1));
vh = P*fh;    % transform   
vh = vh./(lam);   % Solve in Fourier space
vh = P*vh;        % transform back
vh = [0; vh; 0];  % set BC





