Fourier transform for Neumann boundary condition

I need to solve system of two coupled partial differential equations numerically.

$\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}$

The domain of system is a square region.

Boundary condition:

$x = constant \implies \frac{\partial x_1}{\partial x} = \frac{\partial x_2}{\partial x} = 0$

$y = constant \implies \frac{\partial x_1}{\partial y} = \frac{\partial x_2}{\partial y} = 0$

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). Thanks.

-
Boundary conditions prescribed where? (You didn't specify the domain; also, Neumann boundary should have only normal derivative to the boundary being zero; you are specifying the full derivative vanishing...) Also, your question may be better suited for scicomp.stackexchange.com – Willie Wong Sep 10 '12 at 11:11
@WillieWong, Thanks a lot for pointing out scicomp.stackexchange. I will edit question to specify boundary conditions more specifically. – chatur Sep 10 '12 at 15:46
The cosine transform should work. That is, assuming that your square is the unit square $[0,1]\times [0,1]$, expand your function in terms of $\sum_{m,n\in \mathbb{Z}} a_{mn} \cos (\pi m x) \cos(\pi n y)$. – Willie Wong Sep 10 '12 at 16:04