Fourier series coefficients in 2 dimensions

I am trying to solve for the coefficients in the fourier series between $0$ and $\pi$.

$$ \sum_{n=1}^{\infty}\sum_{m=1}^{\infty}A_{nm}\sin(nx)\sin(mx) = 1 $$.

I determined the coefficients to be are as follows.

$$ A_{nm} = \frac{4}{\pi^2}\frac{(1-(-1)^n)(1-(-1)^m)}{nm} $$

However upon putting these into matlab I find that the coefficients are incorrect, as can be seen below.

How exactly would you go about determining the coefficients of this series?

Your expression for the coefficients is fine, up to some constants. Taking $f(x)=1$, we have $$ \int _0^{\pi }dx\int _0^{\pi }dy\, \frac{2 }{\pi }\sin (n x) \sin (m y)f(x) = \frac{2 }{\pi }\frac{(\cos (\pi m)-1) (\cos (\pi n)-1)}{m n} = \frac{2 }{\pi }\frac{1-\left((-1)^m\right) \left(1-(-1)^n\right)}{m n}\,. $$ Then, when we reconstruct a truncation of the Fourier series, i.e., $$ \sum _{n=1}^N \sum _{m=1}^N \frac{2 }{\pi }\frac{1-\left((-1)^m\right) \left(1-(-1)^n\right)}{m n} \frac{2 }{\pi }\sin (n x) \sin (m x)\,, $$ we get an approximation of the constant-1 function. Below is some Mathematica code and the corresponding results. You can see that there is significant Gibbs ringing at the boundaries due to the discontinuity there (due to the mismatch between the boundary conditions of the function and of the basis functions), but otherwise, you can see that it's converging to 1.

This means that you must have something wrong in your MatLab code.

Simplify[Integrate[2/π Sin[n x] Sin[m y], {x, 0, π}, {y, 0, π}], {m, n} ∈ Integers]
u[x_] = Sum[% 2/π Sin[n x] Sin[m y], {n, 1, 20}, {m, 1, 20}];
Plot3D[u[x], {x, 0, π}, {y, 0, π}]