# How to find the variance of a function of a normally distributed random variable. The function is cummulative function of normal distribution

I am wondering how to find the variance of a special function of a normally distributed random variable. More specifically, I am confused by the following question: assume $Y\sim N(\mu,\sigma^2)$, how to caculate $$E(\Phi^2(Y))=\frac{1}{(\sqrt{2\pi})^3\sigma}\int_{-\infty}^{+\infty}\left(\int_{-\infty}^{y}e^{-\frac{z^2}{2}}dz\right)^2 e^{-\frac{(y-\mu)^2}{2\sigma^2}}dy$$ where $\Phi(\cdot)$ is the cumulative function of standard normal distribution? Thank you very much.