Helmholtz FFT solver for 2D fluid simulations in MATLAB

Question

I am performing a 2D fluid simulation in vorticity-streamfunction formulation and instead of a Poisson equation at each time step I need to solve an inhomogeneous Helmholtz equation of the form $$\lambda^2\nabla^2 u = u +f$$ with double-periodic BCs.

I have written the following MATLAB function to do it using FFT

%% Fast solver for 2D inhomogeneous Helmholtz 

function u=helm2d_fft(f,N,h,lambda)

kx = zeros(1, N);
ky = zeros(1, N);
data = zeros(N, N);
data1 = zeros(N, N);
e = zeros(N, N);
u = zeros(N, N);
hh = 2.0*pi/(N*h); % wavenumber indexing

for i = 1:floor(N/2)
    kx(i) = hh*(i-1.0);
    kx(i+floor(N/2)) = hh*(i-floor(N/2)-1);
end
kx(1) = eps;
ky = kx;

for i = 1:N
    for j = 1:N
        data(i,j) = complex(f(i,j), 0.0);
    end
end

e = fft2(data);

e(1,1) = 0.0;
for i = 1:N
    for j = 1:N
        data1(i,j) = e(i,j)/(1-lambda^2*(-kx(i)^2 - ky(j)^2));
    end
end

u = real(ifft2(data1));

end

However it seems that the Helmholtz solver produces a numerical error which is eliminated when I use the following direct solver

% Direct solver

function u=direct_helm_solve(f,N,h,lambda);

F = reshape(f,N*N,1);

A_diag=eye(N)*(1+4*de^2/h^2);
A_diag=A_diag-diag(ones(N-1,1),1)*lambda^2/h^2; 
A_diag=A_diag-diag(ones(N-1,1),-1)*lambda^2/h^2;
A_diag(1,N)=-lambda^2/h^2;
A_diag(N,1)=-lambda^2/h^2;
A_off=-lambda^2*eye(N)/h^2;
A=sparse(N*N,N*N);

for i=1:N
  A((i-1)*N+1:(i-1)*N+N,(i-1)*N+1:(i-1)*N+N)=A_diag; 
end

for i=2:N
  A((i-2)*N+1:(i-2)*N+N,(i-1)*N+1:(i-1)*N+N) =A_off;
 A((i-1)*N+1:(i-1)*N+N,(i-2)*N+1:(i-2)*N+N) =A_off;
end

A(1:N,N*N-N+1:N*N)=-lambda^2/h^2*eye(N);
A(N*N-N+1:N*N,1:N)=-lambda^2/h^2*eye(N);

 uvec = A\F;  
 u=reshape(uvec,N,N);
 end

I can't understand what is wrong with the FFT solver. The Direct solver of course cannot be used for high resolutions.