Evaluating an integral by setting it up as a differential equation

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$ :
• It is written : "the limits being determined by the fact that $I(a)=0$ when $a=0$. I cannot discussed that, considering that it is a specified condition. If it is true, this implies the lower bound of the integral be $0$. The upper bound is necessarily $a$ since we integrate to obtain un function of $a$, not of something else. . – JJacquelin Aug 1 '17 at 16:01
• Right. It can be seen from (2) that $I(a)=0$ for $a=0$. So the only lower bound that works is zero...therefore we have $\int_{0}^{a}\dots$. This is also where the Cauchy PV comes into play since evaluating (2) for $a=0$ requires computing $\mathrm{P.V.\int_{\mathbb{R}}}\dots \mathrm{d}y$ – Aaron Hendrickson Aug 1 '17 at 16:25