Solution
= Solution
Use <tensor-product quadrature> with <composite Simpson's rule> in each coordinate. If there are $m$ subintervals in each direction, then $N=(m+1)^2$ and the <quadrature rule> error is $O(m^{-4})$: apply the one-dimensional <error bound for composite Simpson's rule> successively, using bounded fourth partial <derivatives>. Hence the error is \b[$O(N^{-2})$] in <big O notation>. The rectangular domain makes the product grid convenient.