I was attempting to answer this question, but then I came across a question of my own involving my attempt.
Task: Prove $$\int_0^\infty\frac{\exp(-x^2)}{1+x^2}\mathrm dx=\frac{\pi e}2\text{erfc}(1)$$ Attempt: $$I=\int_0^{\infty}\frac{\exp(-x^2)}{1+x^2}\mathrm dx$$ We then use the Taylor series for the exponential function to find that $$I=\sum_{n\geq0}\frac{(-1)^n}{n!}\int_0^\infty\frac{x^{2n}}{1+x^2}\mathrm dx$$ Setting $x=\tan u$, $$I=\sum_{n\geq0}\frac{(-1)^n}{n!}\int_0^{\pi/2}\tan(u)^{2n}\mathrm{d}u$$ $$I=\sum_{n\geq0}\frac{(-1)^n}{n!}\int_0^{\pi/2}\sin(u)^{2n}\cos(u)^{-2n}\mathrm{d}u$$ And using $$\int_0^{\pi/2}\sin(t)^a\cos(t)^b\mathrm{d}t=\frac{\Gamma(\frac{a+1}{2})\Gamma(\frac{b+1}{2})}{2\Gamma(\frac{a+b}{2}+1)}$$ We have $$I=\frac12\sum_{n\geq0}\frac{(-1)^n}{n!}\Gamma\bigg(\frac{1+2n}{2}\bigg)\Gamma\bigg(\frac{1-2n}{2}\bigg)$$ $$I=\frac12\sum_{n\geq0}\frac{(-1)^n}{n!}\Gamma\bigg(\frac12+n\bigg)\Gamma\bigg(\frac12-n\bigg)$$ Recall the Gamma reflection formula: $$\Gamma(s)\Gamma(1-s)=\pi\csc\pi s\ ,\qquad s\not\in\Bbb Z$$ Since $n\in\Bbb N_0$, we have $\frac12+n\not\in\Bbb Z$, which means we may plug in $s=\frac12+n$: $$I=\frac12\sum_{n\geq0}\frac{(-1)^n}{n!}\pi\csc\bigg(\frac\pi2+\pi n\bigg)$$ $$I=\frac\pi2\sum_{n\geq0}\frac{(-1)^n}{n!}\csc\bigg(\frac\pi2(2n+1)\bigg)$$ Then we recall that $$\sin\bigg(\frac\pi2(2n+1)\bigg)=(-1)^n,\qquad n\in\Bbb Z$$ So we have $$I=\frac\pi2\sum_{n\geq0}\frac{(-1)^n}{n!}\frac1{(-1)^n}$$ $$I=\frac\pi2\sum_{n\geq0}\frac1{n!}$$ $$I=\frac{\pi e}2$$
But $$\frac{\pi e}2\neq \frac{\pi e}2\text{erfc}(1)$$ What did I do wrong? Thanks.
Edit:
I see that $$\int_{\Bbb R^+}\frac{x^{2n}}{1+x^2}\mathrm dx$$ diverges, and as was pointed out in the comments, I can't interchange the $\sum$ and $\int$, but why? The Taylor series converges for all $x\in\Bbb R_0^+$, so what's wrong with the swappage?
Trying to expand the exponential as a series will not work. The integral for each term, except $n=0$, diverges: $$ \sum_{n=0}^\infty\frac{(-1)^n}{n!}\int_0^\infty\frac{x^{2n}}{1+x^2}\,\mathrm{d}x\tag1 $$ The first error in the proposed solution occurs when using $$ \int_0^{\pi/2}\sin(t)^a\cos(t)^b\mathrm{d}t=\frac{\Gamma(\frac{a+1}{2})\Gamma(\frac{b+1}{2})}{2\Gamma(\frac{a+b}{2}+1)}\tag2 $$ This is only valid for $b\gt-1$, so fails for $b=-2n$ when $n\ge1$.
Another approach would be to expand the denominator in a power series. However to address the radius of convergence of the series for $\frac1{1+x^2}$, we need to break the integral into two pieces. $$ \begin{align} \int_0^1\frac{e^{-x^2}}{1+x^2}\,\mathrm{d}x &=\int_0^1\frac{e^{-x^2}}{1+x^2}\,\mathrm{d}x+\int_1^\infty\frac{e^{-x^2}}{1+x^2}\,\mathrm{d}x\tag{3a}\\ &=\sum_{n=0}^\infty(-1)^n\int_0^1e^{-x^2}x^{2n}\,\mathrm{d}x+\sum_{n=0}^\infty(-1)^n\int_1^\infty e^{-x^2}x^{-2n-2}\,\mathrm{d}x\tag{3b} \end{align} $$ As each integral is decreasing and less than $\frac1{2n+1}$, each of these sums converges conditionally. However, I don't see that this can easily be worked into the desired form.
Here is an approach using a trick due to Feynman.
Define $$ F(a)=\int_0^\infty\frac{e^{-ax^2}}{1+x^2}\,\mathrm{d}x\tag4 $$ then $$ F'(a)=-\int_0^\infty\frac{x^2e^{-ax^2}}{1+x^2}\,\mathrm{d}x\tag5 $$ Furthermore, using $D=\frac{\mathrm{d}}{\mathrm{d}a}$, $$ \begin{align} -e^aD(e^{-a}F(a)) &=(I-D)F(a)\tag{6a}\\[8pt] &=\int_0^\infty e^{-ax^2}\,\mathrm{d}x\tag{6b}\\ &=\sqrt{\frac\pi{4a}}\tag{6c} \end{align} $$ Multiplying by $e^{-a}$ and integrating from $b$ to $\infty$ gives $$ \begin{align} e^{-b}F(b)-\overbrace{\lim_{a\to\infty}e^{-a}F(a)}^0 &=\int_b^\infty e^{-a}\sqrt{\frac\pi{4a}}\,\mathrm{d}a\tag{7a}\\ &=\sqrt\pi\int_{\sqrt{b}}^\infty e^{-a^2}\,\mathrm{d}a\tag{7b}\\[3pt] &=\frac\pi2\,\operatorname{erfc}\left(\sqrt{b}\right)\tag{7c} \end{align} $$ Therefore, setting $b=1$, $$ \begin{align} \int_0^\infty\frac{e^{-x^2}}{1+x^2}\,\mathrm{d}x &=F(1)\tag{8a}\\ &=\frac{e\pi}2\,\operatorname{erfc}(1)\tag{8b} \end{align} $$