No.
Using product quadrature rules for multi-dimensional integrals suffers from the so-called curse of dimensionality. An $O(N^{-2})$ accurate rule using N intervals in one-dimension is generally $O(N^{-2})$ accurate when applied as a product rule in $d$ dimensions. However there will be $M = N^d$ evaluations required. So the accuracy is $O(M^{-2/d}).$
Again, if $N$ is the number of intervals in each dimension, then $O(N^{-2})$ accuracy is generally attained assuming the function is well-behaved. Viewed in terms of the total number of quadrature points, however, the order of accuracy gets much worse
To see how the order of accuracy with respect to $N$ carries through consider a simpler mid-point rule.
If $f:[0,1] \rightarrow \mathbb{R}$ and the second-derivative is bounded -- $|f''(x)|\leq K$ -- then the error bound for the mid-point rule applied at the mid-points $\bar{x_i}$ of $N$ uniform invervals is
$$\left|\int_{0}^{1}f(x)\,dx - \frac1{N}\sum_{i=1}^{N}f\left(\bar{x_i}\right)\right|\leq \frac{K}{24N^2}.$$
In the two-dimensional case, $f:[0,1]^2 \rightarrow \mathbb{R}$, where $f$ has bounded second-order partial derivatives we have the same order of accuracy with respect to N. Using a Taylor series approximation, we find that the first- and mixed second-order derivative terms vanish after integrating and the remaining error is
$$\left|\int_{0}^{1}\int_{0}^{1}f(x,y)\,dx\,dy - \frac1{N^2}\sum_{i=1}^{N}\sum_{j=1}^{N}f\left(\bar{x_i},\bar{y_i}\right)\right|\\=\left|\sum_{i=1}^{N}\sum_{j=1}^{N}\int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}\frac1{2}\left[f_{xx}(\xi_x,\xi_y)(x-\bar{x_i})^2+f_{yy}(\eta_x,\eta_y)(y-\bar{y_i})^2\right]\,dx\,dy \right| \\ = N^2O(N^{-4})= O(N^{-2})$$
To provide more detail, using the triangle inequality, the global error bound is
$$E_N \leqslant \sum_{i=1}^{N}\sum_{j=1}^{N}\int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}\frac1{2}\left|f_{xx}(\xi_x,\xi_y)(x-\bar{x_i})^2+f_{yy}(\eta_x,\eta_y)(y-\bar{y_i})^2\right|\,dx\,dy.$$
Assume that the second-order partial derivatives are bounded in absolute value by $K$ and the points are uniformly spaced with $x_i - x_{i-1} = 1/N$ and $y_j - y_{j-1} = 1/N.$
The local error bound on an individual rectangle is
$$E_{ij} = \\ \int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}\frac1{2}\left|f_{xx}(\xi_x,\xi_y)(x-\bar{x_i})^2+f_{yy}(\eta_x,\eta_y)(y-\bar{y_i})^2\right| \, dx \, dy \\ \leqslant \frac{K}{2}\int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}(x-\bar{x_i})^2 \, dx \, dy + \frac{K}{2}\int_{x_{i-1}}^{x_i}\int_{y_{j-1}}^{y_j}(y-\bar{y_i})^2 \, dx \, dy \\ = \frac{K}{2}(y_j -y_{j-1})\int_{x_{i-1}}^{x_i}(x-\bar{x_i})^2 \, dx + \frac{K}{2}(x_i -x_{i-1})\int_{y_{j-1}}^{y_j}(y-\bar{y_i})^2 \, dy \\ = \frac{K}{6N}\left.(x-\bar{x_i})^3\right|_{x_{i-1}}^{x_i} + \frac{K}{6N}\left.(y-\bar{y_j})^3\right|_{y_{j-1}}^{y_j}.$$
With $\bar{x}_i = (x_i + x_{i-1})/2$ and $\bar{y}_i = (y_j + y_{j-1})/2$ we have
$$x_i - \bar{x}_i = -(x_{i-1} - \bar{x}_i) = (x_i - x_{i-1})/2 = 1/2N, \\ y_j - \bar{y}_j = -(y_{j-1} - \bar{y}_j) = (y_j - y_{j-1})/2 = 1/2N,$$
and
$$E_{ij} \leqslant \frac{K}{6N}2\left(\frac{1}{2N}\right)^3 + \frac{K}{6N}2\left(\frac{1}{2N}\right)^3 = \frac{K}{12N^4}.$$
Summing over $N^2$ rectangles, the global error bound is
$$E_N = N^2 \frac{K}{12N^4} = \frac{K}{12N^2}$$
Hence, the local error is $O(N^{-4})$, but the global error is $O(N^{-2})$ as in the one-dimensional case.
Best Answer
You were on the right way: just use the error formula to get
$$h^2 \leq \sqrt{\frac{12 \text{Tol}}{7e}}$$
which yields for $\text{Tol} = 10^{-3}$:
$h \approx 0.025$ and hence $n \approx 40$. All in all, just set $n>40$ and you'll observe an absolute error less that $\text{Tol}$ in your simulation.
The following C++ snippet shows you this fact
If you're not familiar with C/C++, just go to your terminal (in MacOS and Linux it's the same) and type the following in order to compile your source code
Then run it with
to see the output
With $n=40$ I get: