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?). For instance, it has kurtosis of 6 (sharper peak), comparing to 3 for Laplace distribution.

Normal Product Distribution

On the Frequency Function of xy. C.Craig, 1936

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

On the Frequency Function of xy. C.Craig, 1936

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