function phi_galerk_sgpe(phi, V, t, dt, lambda, r, g, k, xstep, T, D)
    noise = sqrt(T * dt / xstep) * sampler_complex(size(phi))
    d = length(size(phi))
    R = CartesianIndices(phi)
    out = similar(phi)
    @. out = (phi - 1im * dt * (g * abs.(phi)^2 + r + lambda * V[1]) * phi + noise)
    out = fft(out)
    for I in R
        out[I] = 1.0 / (1.0 + 1im * dt * sum([k[I[i]]^2 for i = 1:d])) * out[I]
    end
    out = ifft(out)
    @. phi = out
end
