function phi_x = Poisson(rho)
global delx;
eps0 = 10^(-20); %to avoid 1/0


[Nx, Nz] = size(rho); 
rho_hat = fft(rho);
cx = cos((2*pi/Nx)*(0:Nx-1)');
cx = repmat(cx,[1,Nz]);
phi_hat = (rho_hat./(cx-1+eps0))*delx^2/2;
phi = real(ifft(phi_hat));

%error = rho - (circshift(phi,1)+circshift(phi,-1)-2*phi)/delx^2;
%max(sum(error,2)/Nz)

phi_x = (circshift(phi, -1) - circshift(phi,1))/2/delx;

1;


    
