I found a proof that proves \begin{equation} \mathbb{E}[X^{-1}] =\int_{-\infty}^{\infty}\frac{1}{x}\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right)\,\mathrm{d}x =\frac{1}{\sigma^{2}}\exp\left(-\frac{\mu^{2}}{2\sigma^{2}}\right)\int_{0}^{\mu}\exp\left(\frac{t^{2}}{2\sigma^{2}}\right)\,\mathrm{d}t\,, \end{equation} where $X\sim\mathcal{N}(\mu,\sigma^{2})$ and $\mathbb{E}[\cdots]$ is the expectation operator. For those who cannot view the link, I have typeset the proof below.
Question 1: Upon inspection of the integral above (left hand side) it is clear that the integral does not converge. How then can the proof be correct? It would seem that the solution would have to be based on the limit of an improper integral.
Question 2: In the proof below, I follow the author up to Eq. 4 but am not sure how he got the result in Eq. 5. I am not very familiar with differential equations but know enough to observe that the solution in Eq. 5 seems to involve one. Could someone please explain how the result in Eq. 5 is determined?
Proof:
Let, \begin{equation} \mathbb{E}[X^{-1}] =\int_{-\infty}^{\infty}\frac{1}{x}\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right)\,\mathrm{d}x\,. \end{equation}
Substitute $y=x/\sigma$ and introduce $a=\mu/\sigma$ to give \begin{equation} \tag{1} \mathbb{E}[X^{-1}] =\frac{1}{\sigma}\int_{-\infty}^{\infty}\frac{1}{y\sqrt{2\pi}}\exp\left(-\tfrac{1}{2}(y-a)^{2}\right)\,\mathrm{d}y\,. \end{equation}
Define $I(a)$ to be \begin{equation} \tag{2} I(a) =\int_{-\infty}^{\infty}\frac{1}{y\sqrt{2\pi}}\exp\left(-\tfrac{1}{2}(y-a)^{2}\right)\,\mathrm{d}y\,. \end{equation}
Then \begin{equation} \tag{3} \exp\left(\tfrac{1}{2}a^{2}\right)\,I(a) =\int_{-\infty}^{\infty}\frac{1}{y\sqrt{2\pi}}\exp\left(-\tfrac{1}{2}y^{2}+ay\right)\,\mathrm{d}y\,, \end{equation} and \begin{align} \tag{4} \frac{\partial}{\partial a}\left[\exp\left(\tfrac{1}{2}a^{2}\right)\,I(a)\right] &=\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\tfrac{1}{2}y^{2}+ay\right)\,\mathrm{d}y\\ &= \exp\left(\tfrac{1}{2}a^{2}\right)\,. \end{align}
Thus \begin{equation} \tag{5} \exp\left(\tfrac{1}{2}a^{2}\right)\,I(a) =\int_{0}^{a}\exp\left(\tfrac{1}{2}t^{2}\right)\,\mathrm{d}t\,, \end{equation} the limits being determined by the fact that $I(a)=0$ when $a=0$. Therefore \begin{equation} \tag{6} I(a) =\exp\left(-\tfrac{1}{2}a^{2}\right)\int_{0}^{a}\exp\left(\tfrac{1}{2}t^{2}\right)\,\mathrm{d}t\,, \end{equation} and \begin{align} \tag{7} \mathbb{E}[X^{-1}] &=\frac{1}{\sigma^{2}}\exp\left(-\frac{\mu^{2}}{2\sigma^{2}}\right)\int_{0}^{\mu}\exp\left(\frac{t^{2}}{2\sigma^{2}}\right)\,\mathrm{d}t\,,\\ &=\frac{\sqrt{2}}{\sigma}\,\mathcal{D}_{+}\left(\frac{\mu}{\sqrt{2}\sigma}\right)\,. \end{align} where $\mathcal{D}_{+}(z)$ \begin{equation} \mathcal{D}_{+}(z)= e^{-z^{2}}\int _{0}^{z}e^{t^{2}}\,\mathrm{d}t\,, \end{equation} is the Dawson function.
$\int_{-\infty}^{\infty}\frac{1}{x}\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right)\,\mathrm{d}x\quad$ I suppose that you raises the question of convergence for the singular point $x=0$. This would be a serious problem if the lower bound of the integral was $0$. But in the present case, $0$ isn't a bound. The Integral is convergent in the sens of Cauchy Principal Value : http://mathworld.wolfram.com/CauchyPrincipalValue.html
$\tag{4} \frac{\partial}{\partial a}\left[\exp\left(\tfrac{1}{2}a^{2}\right)\,I(a)\right] = \exp\left(\tfrac{1}{2}a^{2}\right)\,.$ Integrating this equation with respect to variable $a$ leads to equation $(5)$, which, on my opinion should be better written with a dummy variable $\alpha$ instead of $a$ :
\begin{equation} \tag{5} \exp\left(\tfrac{1}{2}a^{2}\right)\,I(a) =\int_{0}^{a}\exp\left(\tfrac{1}{2}\alpha^{2}\right)\,\mathrm{d}\alpha\, \end{equation}