$\newcommand{\Cov}{\operatorname{Cov}}$If $X$ and $Y$ are random variables with a bivariate normal distribution and:
May I compute $\Cov(X^m,Y^n)$ for arbitrary positive integers $m$ and $n$?
$\newcommand{\Cov}{\operatorname{Cov}}$Given the new background that the OP provided to his question (most general normal bivariate distribution for variables (X_1, X_2) with non-zero covariance), it is possible to find a general formula for $\Cov(X_1^m, X_2^n)$ as follows:
Consider the generating function $\phi(t_1,t_2)$ for the bivariate distribution given here, in equation (57) (a clean derivation is given for it so I won't repeat it). Then:
$$E(X_1^m X_2^n)=(-i)^{m+n}(\frac{\partial}{\partial t_1})^m(\frac{\partial}{\partial t_2})^n\phi(t_1,t_2)\Big|_{(t_1, t_2)=(0,0)} \tag 1$$
After some tedious and careful algebra we wish to write the generating function $\phi(t_1,t_2)=e^{it_1\mu_1+it_2\mu_2}e^{-\frac{1}{2}(\sigma_1^2t_1^2+2\rho\sigma_1t_1\sigma_2t_2+\sigma_2^2t_2^2)}$ in the form:
$$\phi(t_1,t_2)=e^{-a}e^{-b(t_2+c)^2}e^{-d(t_1+gt_2+h)^2}\tag2$$
which is possible for the values:
$$a=\frac{\mu_1^2}{2\sigma_1^2}+\frac{(\mu_2-\mu_1\rho \frac{\sigma_2}{\sigma_1})^2}{2\sigma^2_2(1-\rho^2)}, \hspace{0.2cm}b=\frac{1}{2}(1-\rho^2)\sigma_{2}^2, \hspace{0.2cm}c=-i\frac{\mu_2-\mu_1\rho \frac{\sigma_2}{\sigma_1}}{2\sigma^2_2(1-\rho^2)}, \hspace{0.2cm}d=\frac{\sigma_1^2}{2},\hspace{0.2cm} g=\frac{\rho\sigma_2}{\sigma_1},\hspace{0.2cm}h=-i\frac{\mu_1}{\sigma_1^2} $$
We change variables $y_2=\sqrt{b}(t_2+c), y_1=\sqrt{d}(t_1+gt_2+h)$ and we find:
$$\begin{align}\frac{\partial}{\partial t_1}&=\sqrt{d}\frac{\partial}{\partial y_1}\\ \frac{\partial}{\partial t_2}&=g\sqrt{d}\frac{\partial}{\partial y_1}+\sqrt{b}\frac{\partial}{\partial y_2}\end{align}$$
Substitute into (1) for the result:
$$E(X_1^mX_2^n)=(-i)^{m+n}(\sqrt{d}\frac{\partial}{\partial y_1})^m(g\sqrt{d}\frac{\partial}{\partial y_1}+\sqrt{b}\frac{\partial}{\partial y_2})^n e^{-a}e^{-y_1^2}e^{-y^2_2}\Big|_{(y_1, y_2)=(c\sqrt{b},h\sqrt{d})}$$
Expanding the parentheses and using the Rodrigues formula for Hermite polynomials we get:
$$E(X_1^mX_2^n)=e^{-a-bc^2-dh^2}i^{m+n}g^n d^{\frac{m+n}{2}}\sum_{k=0}^n\frac{n!}{k!(n-k)!}\Big(\frac{\sqrt{b}}{g\sqrt{d}}\Big)^{n-k}H_{n-k}(h\sqrt{d})H_{m+k}(c\sqrt{b})$$
The expectation values $E(X_1^m)$ and $E(X_2^n)$ can be calculated by setting $n=0$, $m=0$ in the general formula respectively. Also, the imaginary units in the argument of the Hermite polynomials are precisely cancelled out by the imaginary prefactor, so that the final result is real.
EDIT: According to the analysis by @machfour, there is indeed an error in the expression of $c$, which has been corrected for. Due to this, the formula simplifies to the final expression:
$$\small{E(X_1^mX_2^n)=\Big(\frac{i}{\sqrt{2}}\Big)^{m+n}(\rho\sigma_2)^n \sigma_1^m\sum_{k=0}^n {n\choose k} \Big(\sqrt{\frac{1}{\rho^2}-1}\Big)^{n-k}H_{n-k}\Big(\frac{-i\mu_1}{\sigma_1\sqrt{2}}\Big)H_{m+k}\Big(-i\frac{\mu_2\sigma_1-\mu_1\rho}{\sqrt{2}\sigma_1\sigma_2\sqrt{1-\rho^2}}\Big)}$$
I don't have enough reputation to comment, but I think there is a typo/error in DinosaurEgg's answer. I repeated the derivation of (2) in that answer and obtained a slightly different value for $c$ which was off by a factor of 2: \begin{equation*} c = -i \frac{\left(\mu_2 - \mu_1\rho \frac{\sigma_2}{\sigma_1}\right)}{\sigma_2^2(1 - \rho^2)} \end{equation*} while the values I obtained for $a$, $b$, $d$, $g$, and $h$ were the same.
Regardless of the correct value, notice also that since there is no constant term in the exponent of the original characteristic function $$\phi(t_1, t_2) = \exp\left(i t_1 \mu_1 + i t_2 \mu_2 - \frac{1}{2}(t_1^2 \sigma_1^2 + 2t_1 t_2 \sigma_1 \sigma_2 \rho + t_2^2 \sigma_2^2) \right),$$ examining (2) suggests that we should have $$ - a - bc^2 - dh^2 = 0$$ as an identity. Thus the final formula simplifies to \begin{align*} E(X_1^mX_2^n) &=i^{m+n}g^n d^{\frac{m+n}{2}}\sum_{k=0}^n\frac{n!}{k!(n-k)!}\Big(\frac{\sqrt{b}}{g\sqrt{d}}\Big)^{n-k}H_{n-k}(h\sqrt{d})H_{m+k}(c\sqrt{b}) \\ &=\left(\frac{i}{\sqrt{2}}\right)^{m+n} \sigma_1^m (\rho \sigma_2)^n \sum_{k=0}^n\frac{n!}{k!(n-k)!}\Big(\sqrt{1/\rho^2 - 1}\Big)^{n-k}H_{n-k}(-\frac{i \mu_1}{\sqrt{2} \sigma_1})H_{m+k}(c\sqrt{b}) \\ \end{align*}
You have said$\newcommand{\Cov}{\operatorname{Cov}}$ \begin{align} & X\sim N(\mu_X,\sigma^2_X), \\[4pt] & Y \sim N(\mu_Y,\sigma^2_Y), \\[4pt] & \Cov(X,Y) = 0. \end{align} That falls short of specifying the joint distribution of $(X,Y).$ If it were further specified that $X,Y$ are jointly normally distributed, then the covariance can be $0$ only if $X,Y$ are independent. If $X,Y$ are independent, then so are $X^m,Y^n,$ so their covariance is also $0.$
Here is a simple example: Suppose $X\sim N(\mu_X, \sigma_X^2),$ let $Z = \text{the “z-score''} = (X-\mu_X)/\sigma_X,$ and independently of $X$ you toss a coin. Then let $Y = \mu_Y \pm \sigma_Y Z, $ where the choice between $\text{“}\pm\text{''}$ is determined by the coin toss. Then $X,Y$ have covariance $0$ and have just the distributions you specified in the question, but they are NOT JOINTLY normally distributed and not independent.
But if you assume joint normality, which in the case of two random variables means bivariate normality, then the answer is just as in the first paragraph above.
