Matlab's fft2 is a wrapper around fftn, but the order of operation of 1D fft along dimensions differs. To cut to the chase, the following gives a non-zero answer. Why is this? I can find no bug reports online.
foo=rand(64);
err=sum(sum(fft(fft(foo,[],2),[],1)-fft(fft(foo,[],1),[],2)))