The following question is a simple case of a type of problem that occurs in combinatorial enumeration problems.
Define $$F(x_1,\ldots,x_n) = \frac{1}{(2\pi)^{n/2}}\exp\biggl( -\frac12\sum_{j=1}^n x_j^2 + i\,n^{-1/4} \sum_{j=1}^n x_j^3 \biggr), $$ where $i$ is the imaginary unit. Let $\epsilon\gt 0$ be small enough. Let $$I_0(n) = \int_{-n^\epsilon}^{n^\epsilon}\cdots \int_{-n^\epsilon}^{n^\epsilon} F(x_1,\ldots,x_n) ~ dx_1\cdots dx_n. $$ We can estimate $I_0(n)$ as $n\to\infty$ by factoring $$F(x_1,\ldots,x_n) = \prod_{j=1}^n \frac{1}{\sqrt{2\pi}} \exp\bigl(-x_j^2/2+in^{-1/4}x_j^3\bigr),$$ which separates the integral into a product of $n$ 1-dimensional integrals. This very easily gives $$I_0(n) = \exp\bigl( -n^{1/2}/2 + O(n^{-3})\bigr).$$ So far so good. But now let $A=A(n)$ be a symmetric positive definite matrix. We can assume it is pretty nice, say all eigenvalues bounded between two positive constants independent of $n$. Now define $$G(x_1,\ldots,x_n) = \frac{1}{(2\pi)^{n/2}}\exp\biggl( -\frac12 \mathbf{x}^TA\mathbf{x} + i\,n^{-1/4} \sum_{j=1}^n x_j^3 \biggr). $$ How do we estimate $$I_1(n) = \int_{-n^\epsilon}^{n^\epsilon}\cdots \int_{-n^\epsilon}^{n^\epsilon} G(x_1,\ldots,x_n) ~ dx_1\cdots dx_n? $$ The source of the difficulty is that $\int G$ is exponentially smaller than $\int |G|$, so as soon as you approximate the integrand the answer goes away. What to do?
Note that in this case (and very commonly in practice) symmetry implies that $I_1(n)$ is real. So we can discard the imaginary part of the integrand and integrate only the real part, which is $$ \mathfrak{R}G(x_1,\ldots,x_n) = \frac{1}{(2\pi)^{n/2}}\exp\biggl( -\frac12 \mathbf{x}^TA\mathbf{x}\biggr)\cos\biggl(n^{-1/4} \sum_{j=1}^n x_j^3 \biggr). $$ Not sure that helps.
ADDED 2 Jan 2018. Mikhail Isaev and I now have a method that can precisely estimate such integrals in many general cases. It isn't completely written up yet.