Arkadiusz gives the answer in the case of two independent Gaussians. A simple technique to reduce the correlated case to the uncorrelated is to diagonalize the system. The intuition which I use is that for two random variables, we need two "independent streams of randomness," which we then mix to get the right correlation structure.

Let $X \sim N(0,\sigma_x)$ and let $Z \sim N(0,1)$ be two independent normals. Define

$Y = \tfrac{\rho \sigma_y}{\sigma_x} X + \sqrt{1-\rho^2}\sigma_y Z$.

Check that $\mathbb E Y^2 = \sigma_y^2$ and $\mathbb E XY = \rho \sigma_x \sigma_y$; this completely determines the bivariate Gaussian case you're interested in.

Now, $XY = \tfrac{\rho \sigma_y}{\sigma_x} X^2 + \sqrt{1-\rho^2}\sigma_y XZ$. The $X^2$ part has a $\chi^2$-distribution, familiar to statistics students; the $XZ$ part is comprised of two independent Gaussians, hence Arkadiusz's answer gives the distribution of that random variable. Let $f_{X^2}$ and $f_{XZ}$ be the two density functions, respectively. The density of a sum of two random variables is the convolution of their density functions: $f_{XY} = f_{X^2} * f_{XZ}$.

Arkadiusz gives the answer in the case of two independent Gaussians. A simple technique to reduce the correlated case to the uncorrelated is to diagonalize the system. The intuition which I use is that for two random variables, we need two "independent streams of randomness," which we then mix to get the right correlation structure.

Let $X \sim N(0,\sigma_x)$ and let $Z \sim N(0,1)$ be two independent normals. Define

$Y = \tfrac{\rho \sigma_y}{\sigma_x}X + \sqrt{1-\rho^2}\sigma_y Z$.

Check that $\mathbb E Y^2 = \sigma_y^2$ and $\mathbb E XY = \rho \sigma_x \sigma_y$; this completely determines the bivariate Gaussian case you're interested in.

Now, $XY = \tfrac{\rho \sigma_y}{\sigma_x} X^2 + \sqrt{1-\rho^2}\sigma_y XZ$. The $X^2$ part has a $\chi^2$-distribution, familiar to statistics students; the $XZ$ part is comprised of two independent Gaussians, hence Arkadiusz's answer gives the distribution of that random variable.

Edit: As Robert Israel points out in the comments, I made a mistake in my final conclusion: the random variables $X^2$ and $XZ$ are uncorrelated, though certainly not independent. Nonetheless, the problem is essentially resolved at this point, since we have reduced the problem of understanding the product $XY$ to a sum of uncorrelated random variables $X^2$ and $XZ$ with known distributions.