It seems that no closed-form expression is known for the correlated case.

For the uncorrelated case $\rho=0$ the distribution of XY is $\frac{1}{\pi \sigma_x \sigma_y}K_0(\frac{|u|}{\sigma_x \sigma_y})$, where $K_0(x)$ is the modified Bessel function of the second kind. This distribution differs from the distribution you gave (is it Laplace distribution?). A notable difference: it has kurtosis of 6 (sharper peak), comparing to 3 for Laplace distribution.

[Normal Product Distribution][1]

[On the Frequency Function of xy. C.Craig, 1936][2]


  [1]: http://mathworld.wolfram.com/NormalProductDistribution.html
  [2]: http://web.ipac.caltech.edu/staff/fmasci/home/statistics_refs/GaussianProduct1934.pdf

Edit: answering the second question, the method is valid, but my guess is that the Jacobian determinant was skipped while doing the substitution.