Why am I not getting the intended coefficients in this 2D Fourier demonstration?
Your misconception is to assume that you can compute Fourier coefficients with the fft. You cannot (really). Fourier coefficients are found by integrating a continuous time (periodic) function after multiplying it with a harmonic kernel. The cosine function gives two real-valued Fourier coefficients equal to half its amplitude in this case. 

You are considering a matrix though which means that you have sampled your cosine function at a discrete grid (which is not a problem in your case) and truncated it to a finite width (which is exactly the problem you have). Then you are computing the FFT which is a fast way of computing a DFT (Discrete Fourier Transform). The DFT (of discrete vectors) is not the same as a Fourier transform (which is defined on continuous functions).

That said, your DFT can approximate Fourier coefficients, as the Fourier integral can be approximated by a sum, which is precisely the definition of the DFT then. The approximation errors that happen in this case are sampling and truncation.

For functions that are strictly bandlimited (like your cosine), sampling does not cause any error as long as you satisfy the Nyquist criterion. Conversely, for functions that are strictly time-limited, truncation does not cause any error as long as the whole function fits in your window. Unfortunately, it can be shown that no function can be time-limited and band-limited at the same time. For this reason, the DFT can never accurately predict the result of the Fourier transform.

There is hope though, since you're not after a Fourier transform but instead you're after Fourier coefficients. These are defined on periodic continuous-time functions. For these, the DFT can give you accurate predictions of the Fourier coefficients, as long as you make sure that the samples you are considering contain a whole period of the signal. This needs to be accurate, i.e., if you periodify your extracted window it should look like the original function. For instance, if the period of your function is 1 and you want to extract four samples, they need to be at 0, 0.25, 0.5, 0.75 (i.e., not 0, 1/3, 2/3, 1, which is a very common mistake in this setting).

Your Fourier coefficients will then be real-valued if and only if your function is symmetric, which in the sampled case means x[n] = x[N-n] for n=1, 2, ..., N-1.  If you sample your cosines correctly, this should work.