function [u,x]=FDPoisson1d(f,a,b,n);
% FDPOISSON1D solves the 1d Poisson equation with periodic boundary conditions
%   u=Poisson1d(f,a,b,n); solves the 1d poisson equation with right
%   hand side function f and periodic boundary conditions using centered finite
%   differences and n+1 mesh points. Note that f is made to have mean zero.

h=(b-a)/n; x=(a:h:b)';                        % finite difference mesh
e=ones(n,1);
A=spdiags([e -2*e e],-1:1,n,n)/h^2;           % finite difference matrix
A(n,1)=1/h^2;
b=feval(f,x(1:n));                            % check mean in the right
bs=sum(b);                                    % hand side and remove if
if abs(bs)>100*eps,                           % necessary
  warning(['right hand side function mean ' num2str(bs) 'removed']);
  b=b-bs/n;
end;
b(1)=0;                                       % replace redundant equation
A(1,1:n)=ones(1,n);                           % by mean zero constraint
u=A\b;
u(n+1)=u(1);                                  % solution is periodic
