Following @whuber suggestion, using the [quotient][1] formula:

$X:=Z/U$

The PDF: $f_X(x)=\int_{u=0}^1 1\cdot \frac{1}{\sqrt{2\pi}}e^{-(xu)^2/2}du = \frac{1}{\sqrt{2\pi}}[-\frac{e^{-x^2u^2/2}}{x^2}]_{u=0}^1 = \frac{1}{x^2}(\phi(0)-\phi(x))$

The CDF: $\int_{x=-\infty}^tf_X(x)dx= \frac{1}{\sqrt{2\pi}}\int_{x=-\infty}^t(\frac{1}{x^2}-\frac{1}{x^2}e^{-x^2/2})dx=^* \frac{1}{\sqrt{2\pi}}(-\frac{1}{t}+\frac{1}{t}e^{-t^2/2}+\int_{-\infty}^t e^{-t^2/2}dt )$

$= \Phi(t)-\frac{\phi(0)-\phi(t)}{t}$

Where I used integration by parts to calculate the $*$ part.


Another way is directly through the CDF:

$F_X(x)=P(X\le x) = P(Z/U \le x) = \int_{u=0}^1 P(Z\le xu)du = \int_{u=0}^1 \frac{1}{2}(1+erf(xu/\sqrt2))du = 0.5 + 0.5\int_{u=0}^1erf(xu/\sqrt2)du$

Replace $xu/\sqrt2 = t \Rightarrow du = \frac{\sqrt2}{x} dt$

$\int_{u=0}^1erf(xu/\sqrt2)du = \frac{\sqrt2}{x}\int_{t=0}^{\frac{x}{\sqrt2}}erf(t)dt = \frac{\sqrt2}{x}[t\cdot erf(t) + \frac{e^{-t^2}}{\sqrt \pi}]_0^{\frac{x}{\sqrt2}}=$

$erf(x/\sqrt2) + \frac{\sqrt2}{x}\frac{e^{-x^2/2}}{\sqrt \pi}-\frac{\sqrt2}{x\sqrt\pi}$

$F_X(x)=0.5 + 0.5erf(x/\sqrt2)+\frac{\phi(x)-\phi(0)}{x} = \Phi(x)+\frac{\phi(x)-\phi(0)}{x} $

And then the PDF can be obtained from taking the derivative.


  [1]: https://www.jstor.org/stable/2235953?seq=1