The PDF of the random variable Z=Y+Y*X I have two independent random variables, $Y$ and $X,$ where
$Y$ is a random variable with a Gaussian distribution and a zero mean. 
$X$ is a random variable with a Gaussian distribution and a zero mean. 
I then create a new random variable:
$$Z = Y+Y*X$$
Is there any way to determine the PDF of $Z$ and how can I compute the correlation between $Y$ and $YX$ if it exists?

Let $U = (X+Y)/\sqrt{2}$ and $V = (X-Y)/\sqrt{2}.$  Because these are linear combinations of the bivariate Normal variable $(X,Y),$ $(U,V)$ has a bivariate Normal distribution.  Direct calculation of the means and covariances shows $U$ and $V$ are uncorrelated standard Normal variables, whence they are independent.
Algebra shows that
$$Z = Y + YX = ((U + \sqrt{1/2})^2 - (V + \sqrt{1/2})^2)/2,$$
a difference of two independent identically distributed terms.  Each of these terms has--by definition--a noncentral $\chi^2(1,1/2)$ distribution..  Consequently, writing $\psi$ for the common characteristic function, the c.f. of $Z$ must be
$$\psi_{Z}(t) = \psi(t/2)\psi(-t/2) = \frac{\exp\left(\frac{it/4}{1-it}\right)}{\sqrt{1-it}} \frac{\exp\left(\frac{-it/4}{1+it}\right)}{\sqrt{1+it}} = \frac{1}{\sqrt{1 + t^2}} \exp\left(\frac{-t^2/2}{1+t^2}\right).$$
The density of $Z$ is the inverse Fourier transform of $\psi_{Z}.$  It can be computed as a cosine transform because $Z$ is (obviously) symmetric,
$$f_{Z}(z) = \frac{1}{2\pi} \int_{0}^{\infty} \frac{2 \cos(t z)}{\sqrt{1 + t^2}} \exp\left(\frac{-t^2/2}{1+t^2}\right)\,\mathrm{d}t.$$
100,000 simulated values of $X$ and $Y$ produced this histogram of $Z.$  On it I have plotted, in red, the values of $f_Z$ computed by numerically integrating this formula.

The integration could be a little dicey because the integral, as a Riemann integral, does not converge absolutely.  However, it does converge and, as it turns out, the calculation is not onerous.  A PC workstation can compute several thousand values a second.
Because even substantial algebraic errors in the calculation could still produce a close agreement in this plot, I simulated a million values of $Z,$ of which this next histogram is a detail near the origin where the integration becomes most difficult:

The agreement remains excellent.

BTW, because $X$ and $Y$ have finite variance, the covariance of $Y$ and $YX$ exists and is finite.  Moreover, since $X$ and $-X$ have the same distribution,
$$\operatorname{Cov}(Y, YX) = \operatorname{Cov}(Y, Y(-X)) = -\operatorname{Cov}(Y, YX)$$
implies this covariance must be $0.$
