function ft = f(at)
    global kvec c epsilon
    filt=[ones(1,32),zeros(1,65),ones(1,31)];
    u=fft(at.*filt);
    u2=u.*u/2;
    b=ifft(u2);
    ft=-kvec.*kvec*epsilon.*at+1i*kvec.*b;
end

