FFTN on 2D distribution

I got a known 2D distrubution of a magnetic field and and the following equations to expand my field with 3D Fourier Series:

and are the amplitudes of the harmonics, and are the phase angles. After some research I found matlab's fftn command. But I dont quite understand how I need to prepare my data and so I can feed it properly to the algorithm and get results for and .

Let's assume I got a 2D-matrix with amplitudes of my magnetic field strength ( and ):

S = rand(10,10)

S = 10×10

0.0774 0.0586 0.1233 0.8253 0.4141 0.3026 0.0886 0.1782 0.8257 0.9804 0.8358 0.5366 0.9505 0.7182 0.9771 0.7725 0.3247 0.0484 0.0352 0.1475 0.3372 0.9373 0.7688 0.6740 0.9309 0.3293 0.9938 0.3643 0.1486 0.8021 0.9284 0.3114 0.5688 0.9209 0.5394 0.2210 0.3374 0.2097 0.7423 0.9325 0.6439 0.4171 0.2081 0.5371 0.6221 0.8705 0.3285 0.2187 0.2559 0.0886 0.5821 0.6240 0.9384 0.6517 0.3809 0.6432 0.6882 0.4982 0.1777 0.0094 0.7878 0.7093 0.8652 0.9309 0.0106 0.4748 0.7695 0.4599 0.5664 0.9333 0.0223 0.6135 0.9721 0.2218 0.9598 0.0652 0.8725 0.6674 0.9581 0.1444 0.9269 0.1790 0.5003 0.2586 0.5871 0.8213 0.7123 0.4241 0.7743 0.1457 0.4185 0.4395 0.3895 0.7202 0.0797 0.7705 0.2463 0.8244 0.3332 0.1039

Since The output Y is the same size as X according to the description, I need to modifiy the matrix first. Is this done via a meshgrid or do I understand something wrong?