I'd like to ask about time evolution and FFTW 3d:
I have prepared initial gaussian packet : psi(x,y,z) = exp(-(x^2+y^2+z^2)).
I'm interested in numerical time evolution, so I prepared 3d forward FFT using FFTW, then multiply psi(k_x,k_y,k_z) by coeffitient and do backward fft. But unfortuanetly something goes wrong and I don't know why.
There is part of this code:

allocate(psi_3d_x(N_x,N_y,N_fi),psi_3d_k(N_x,N_y,N _fi))
w_x = 2*pi**2/DBLE(M)/dx**2/DBLE(N_x)**2.d0*dt
w_y = 2*pi**2/DBLE(M)/dy**2/DBLE(N_y)**2.d0*dt
w_fi = 2*pi**2/DBLE(M)/dfi**2/DBLE(N_fi)**2.d0*dt
call dfftw_plan_dft_3d(plan_3d_f,N_x,N_y,N_fi,psi_3d_x, psi_3d_k,FFTW_FORWARD,FFTW_estimate)
call dfftw_plan_dft_3d(plan_3d_b,N_x,N_y,N_fi,psi_3d_k, psi_3d_x,FFTW_backward,FFTW_estimate)

do i = 1, N_x
do j = 1, N_y
do d = 1, N_fi
x = -L/2.d0 + (i-1)*dx
y = -L/2.d0 + (j-1)*dy
fi = -L/2.d0 + (d-1)*dfi
psi_3d_x(i,j,d) = dcmplx(exp(-DBLE(M)*omega*(abs(x)**2+abs(y)**2+abs(fi)**2)/2.d0)*(DBLE(M)*omega/pi)**0.75d0,0.d0)
end do
end do
end do

call dfftw_execute_dft(plan_3d_f,psi_3d_x,psi_3d_k)


do i = 2,N_x/2
do j = 2,N_y/2
do d = 2, N_fi/2
w_x_1 = w_x*DBLE(i-1)**2
w_y_1 = w_y*DBLE(j-1)**2
w_fi_1 = w_fi*DBLE(d-1)**2
psi_3d_k(i,j,d) = psi_3d_k(i,j,d)*dcmplx(cos(w_x_1),-sin(w_x_1))*dcmplx(cos(w_y_1),-sin(w_y_1))*dcmplx(cos(w_fi_1),-sin(w_fi_1))
end do
end do
end do
do i = N_x/2+1,N_x
do j = N_y/2+1,N_y
do d = N_fi/2+1,N_fi
w_x_1 = w_x*DBLE(N_x-(i-1))**2
w_y_1 = w_y*DBLE(N_y-(j-1))**2
w_fi_1 = w_fi*DBLE(N_fi-(d-1))**2
psi_3d_k(i,j,d) = psi_3d_k(i,j,d)*dcmplx(cos(w_x_1),-sin(w_x_1))*dcmplx(cos(w_y_1),-sin(w_y_1))*dcmplx(cos(w_fi_1),-sin(w_fi_1))
end do
end do
end do

call dfftw_execute_dft(plan_3d_b,psi_3d_k,psi_3d_x).


Can anyone help me ?
Thanks a lot,
Marcin