I need to solve system of two coupled partial differential equations numerically. \begin{align}\frac{\partial x_1}{\partial t} &= c_1\nabla ^2 x_1 + f_1(x_1,x_2)\\ \frac{\partial x_2}{\partial t} &= c_2\nabla ^2 x_2 + K\frac{\partial x_1}{\partial t}\end{align}The domain of system is a square region. Boundary condition:\begin{align}x &= \text{constant} \implies \frac{\partial x_1}{\partial x} = \frac{\partial x_2}{\partial x} = 0\\y &= \text{constant} \implies \frac{\partial x_1}{\partial y} = \frac{\partial x_2}{\partial y} = 0\end{align}I tried to solve this system with Fourier transform. Solution becomes unstable after few iterations. I have solved this system earlier with finite difference scheme and it worked well so I know that constants of system are perfectly fine.My question is can Fourier transform be used to solve these equations? I read somewhere that it because of Neumann boundary condition one cannot apply Fourier transform.  Is this correct? If yes what is alternative?(I have read that cosine transform should be used but want to confirm). 
Fourier transform for Neumann boundary condition
numerics;pde;fourier analysis