Let $X$ and $Y$ be independent random variables with nonzero variances. I'm looking to find the correlation coefficient $\rho$ of $Z=XY$ and $X$ in terms of the means and variances of $X$ and $Y$, i.e. $\mu_X, \mu_Y, \sigma^2_X, \sigma^2_Y$.
(I have searched different methods online, including Correlation between X and XY. However, I'm wondering if I could use a simple calculation-approach rather than using moments as well.)
The result I obtained, along with the steps I've used, is the following:
$$ \begin{align} \rho & = \frac{\text{Cov}(Z,X)}{\sigma_Z\sigma_X}\\[1em] & = \frac{E\left[\left(Z-\mu_Z\right)\left(X-\mu_X\right)\right]}{\sigma_Z\sigma_X} \\[1em] & = \frac{E\left[\left(XY-\mu_X\mu_Y\right)\left(X-\mu_X\right)\right]}{\sqrt{E\left[\left(XY\right)^2\right]-\left[E\left(XY\right)\right]^2}\cdot\sigma_X} \\[1em] & = \frac{E\left(X^2Y\right)-\mu_X^2\mu_Y}{\sqrt{E\left(X^2\right)E\left(Y^2\right)-\left[E\left(X\right)\right]^2\left[E\left(Y\right)\right]^2}\cdot\sigma_X} \\[1em] & = \frac{E\left(X^2\right)E\left(Y\right)-\mu_X^2\mu_Y}{\sqrt{\left(\sigma_X^2+\mu_X^2\right)\left(\sigma_Y^2+\mu_Y^2\right)-\mu^2_X\mu^2_Y}\cdot\sigma_X} \\[1em] & = \frac{\mu_Y\left[E\left(X^2\right)-\mu^2_X\right]}{\sqrt{\sigma^2_X\sigma^2_Y+\sigma_X^2\mu_Y^2+\sigma_Y^2\mu_X^2}\cdot\sigma_X} \\[1em] & = \frac{\mu_Y\sigma_X^2}{\sqrt{\sigma^2_X\sigma^2_Y+\sigma_X^2\mu_Y^2+\sigma_Y^2\mu_X^2}\cdot\sigma_X} \\[1em] & = \frac{\mu_Y\sigma_X}{\sqrt{\sigma^2_X\sigma^2_Y+\sigma_X^2\mu_Y^2+\sigma_Y^2\mu_X^2}} \end{align} $$
which is seemingly different from the result from the moment approach used in Correlation between X and XY. In which step has an error in my calculation occurred (if any), and how can I obtain $\rho$ from the approach I am trying to use?