Following @whuber suggestion, using the quotient 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.
I think the CDF can also be understood geometrically, looking at the joint PDF, for a positive $X=x$ we are looking for the red area (the PDF is the 3rd dimension above this area):
We can calculate this as the area of the rectangle at the point where $X=z, u=1$, which is just $\Phi(x)$, minus the area of the green triangle, which is proportional to $\phi(0)-\phi(x)$ by $x$, since the bigger the $x$ value, the more you replace high density for low density.
i.e. when you increase $x$ you give up the yellow triangle that is high density, for the cyan triangle and rectangle which isare low density.
To see this algebraically, we are interested in the volume under the green triangle area:
$\int_{z=0}^x\int_{u=0}^{z/x}\phi(z)dudz = \int_{z=0}^x\phi(z)\frac{z}{x}dz=\frac{1}{x}(-\phi(x)+\phi(0))$