Let's say the two variables are x and y, and we want an FFT of f(x,y).
1) Build a NxM array of samples of f(x,y), N samples one dimension, M in the other.
2) Do N FFTs (one on each row), each M long, save the results back into the same array.
3) Do M FFTs (one on each column), each N long on the result of step 2.
4) That's it, you are done! Surprisingly simple, huh?
In Sonnet, we use a 2-D FFT to calculate the coupling between subsections. The fields are represented as a sum of waveguide modes. Waveguide modes have mode numbers m, and n (like TM 01, etc.). The coupling between subsections (where current on one subsection generates voltage on another) is a 2-D sum of the sines and cosines of the waveguide modes. Since we are just summing the sines and cosines (not complex exponentials, as is done in the FFT), we only have to do about half the work of a complete complex FFT.
The FFT is why the Sonnet EM analysis has such a large dynamic range. The coupling between subsections is calculated with an FFT to full numerical precision, no numerical integration error. It is also why we do so well when there are a larger number of layers, just add on a couple more FFTs (which typically only require a couple more seconds in most cases).