Let $f(x,y)$ be a Gaussian pdf for some known mean and covariance.
Given $(x_0, y_0)$ and $(x_N, y_M)$ such that $$\int_{x_0}^{x_N} \int_{y_0}^{y_M} f(x,y) dy dx \approx 1$$
I would like to split the intervals $[x_0, x_N]$ and $[y_0, y_M]$ into $N+1$ and $M+1$ respectively (i.e. $\{x_0, x_1, \dots x_N\}, \{y_0, y_1, \dots y_M\}$) such that $$\int_{x_i}^{x_{i+1}} \int_{y_{j}}^{y_{j+1}} f(x,y) dy dx \approx \frac{1}{NM}$$
Work So Far My initial approach was to assume we had $N,M$ large enough to use the approximation
$$\int_{x_i}^{x_{i+1}} \int_{y_{j}}^{y_{j+1}} f(x,y) dy dx \approx f(x_i,y_j)(x_{i+1}-x_i)(y_{j+1}-y_j)$$
We then want to solve the system of equations $$f(x_i,y_j)(x_{i+1}-x_i)(y_{j+1}-y_j) = \frac{1}{NM}$$
This gives $N+M-4$ unknowns and $NM$ equations. I was thinking of some sort of iterative method to solve this but was not sure of the best way to formulate it. Note that in the 1D case we can just substitute in our $x_0$ value and loop over all $i$. How best to proceed?
Start from the center and divide each semi interval in $N/2$ (with $N$ even).
so $x_0 = \mu_x$ and the interval go from $x_{-N/2}...x_0...x_{N/2}$
You can then use the inverse error function to compute each $x_i$ such that
$$\frac{1}N=\frac{1}{2}\left(erf\left(\frac{x_{i+1}-\mu_x}{\sqrt{2\sigma_x^2}}\right)-erf\left(\frac{x_{i}-\mu_x}{\sqrt{2\sigma_x^2}}\right)\right)$$
similarly over the y axis.