Mathematica gives the result of the above integral in terms of a Meijer G function, or numerically
$$ PV\int_{0}^{\infty} dx \frac{e^{-x^2}}{(x-1)^2} \approx -0.168 - 2.31 i $$
Ideally I would like to be able to approximate this answer using some numerical integration. But I simply can't understand how there is a finite answer for this. Around the singularity we have
$$ \frac{e^{-x^2}}{(x-1)^2}\Bigg|_{x=1-\epsilon} \approx \frac{1}{e}\Big(\frac{1}{\epsilon^2} - \frac{2}{\epsilon}+1+O(\epsilon) \Big)$$ $$ \frac{e^{-x^2}}{(x-1)^2}\Bigg|_{x=1+\epsilon} \approx \frac{1}{e}\Big(\frac{1}{\epsilon^2} + \frac{2}{\epsilon}+1+O(\epsilon) \Big)$$
I can understand how the first order terms will 'cancel' on the left and right sides of the singularity in a principal value sense, but the leading terms both go to infinity in the same direction.
I would appreciate any guidance with understanding this. Thanks!
Edit: The exact expression from Mathematica is
-(MeijerG[{{0, 1/2, 1}, {}}, {{1/2, 1}, {}}, -1, 1/2]/\[Pi])
which I believe is $$ -\frac{1}{\pi} G_{3,2}^{2,3}\Big(-1,\frac{1}{2}\Big| \begin{matrix}0,1/2,1\\1/2,1\end{matrix} \Big) $$
I was doing the same when @Felix Marin's comment appeared and the result is $$\frac{\pi \, \text{erfi}(1)+\text{Ei}(1)-e \left(1+\sqrt{\pi }\right)}{e}\sim-0.167837 $$