Detailing Srivatsan Narayanan's solution. It is known that the functional equation of the gamma function may be derived applying the integration by parts technique. Its value at $1/2$ may be evaluated by computing a double integral over the first quadrant in Cartesian and polar coordinates. Let's apply similar ideas in this case. Let $f(x)=x^{2}e^{-x^{2}}$. Since $ f(-x)=f(x)$ the integral $\int_{-\infty }^{\infty }f(x)\mathrm{d}x=2\int_{0}^{\infty }f(x)\mathrm{d}x$. Integrating by parts $ \int x^{2}e^{-x^{2}}\mathrm{d}x=\int x\cdot xe^{-x^{2}}\mathrm{d}x, $ since $ \int xe^{-x^{2}}\mathrm{d}x=-\frac{1}{2}e^{-x^{2}} $ and $\frac{dx}{dx}=1$, we get $\begin{eqnarray*} \int x^{2}e^{-x^{2}}\mathrm{d}x &=&x\left( -\frac{1}{2}e^{-x^{2}}\right) -\int -\frac{ 1}{2}e^{-x^{2}}\,\mathrm{d}x \\ &=&-\frac{1}{2}xe^{-x^{2}}+\frac{1}{2}\int e^{-x^{2}}\,\mathrm{d}x.\tag{0} \end{eqnarray*}$
And so, $\begin{eqnarray*} \int_{0}^{\infty }x^{2}e^{-x^{2}}dx &=&\left. -\frac{1}{2} xe^{-x^{2}}\right\vert _{0}^{\infty }+\frac{1}{2}\int_{0}^{\infty }e^{-x^{2}}\,\mathrm{d}x \\ &=&\left( \lim_{c\rightarrow \infty }-\frac{1}{2}ce^{-c^{2}}\right) +\frac{1 }{2}0e^{-0^{2}}+\frac{1}{2}\int_{0}^{\infty }e^{-x^{2}}\,\mathrm{d}x \\ &=&0+0+\frac{1}{2}\int_{0}^{\infty }e^{-x^{2}}\,\mathrm{d}x \\ &=&\frac{1}{2}\int_{0}^{\infty }e^{-x^{2}}\,\mathrm{d}x.\tag{1} \end{eqnarray*}$
Consequently, $ I:=\int_{-\infty }^{\infty }x^{2}e^{-x^{2}}\mathrm{d}x=2\int_{0}^{\infty }x^{2}e^{-x^{2}}\mathrm{d}x=\int_{0}^{\infty }e^{-x^{2}}\,\mathrm{d}x.\tag{2} $
To evaluate this last integral we compute the following double integral in Cartesian and polar coordinates ($r^{2}=x^{2}+y^{2}$, $x=r\cos \theta ,y=r\sin \theta $). Since the Jacobian of the transformation $\frac{\partial (x,y)}{\partial (r,\theta )}=r$, we have $ \int_{x=0}^{\infty }\int_{y=0}^{\infty }e^{-x^{2}-y^{2}}\mathrm{d}x\mathrm{d}y=\int_{\theta =0}^{\pi /2}\int_{r=0}^{\infty }e^{-r^{2}}r\mathrm{d}r\mathrm{d}\theta.\tag{3} $ Comparing the LHS $ \int_{x=0}^{\infty }\int_{y=0}^{\infty }e^{-x^{2}-y^{2}}\mathrm{d}x\mathrm{d}y=\left( \int_{0}^{\infty }e^{-x^{2}}\mathrm{d}x\right) \left( \int_{0}^{\infty }e^{-y^{2}}\mathrm{d}y\right) =I^{2}\tag{4} $ with the RHS $\begin{eqnarray*} I^2 &=&\int_{\theta =0}^{\pi /2}\int_{r=0}^{\infty }e^{-r^{2}}r\mathrm{d}r\mathrm{d}\theta = \frac{\pi }{2}\int_{0}^{\infty }e^{-r^{2}}r\mathrm{d}r \\ &=&\frac{\pi }{2}\left. \left( -\frac{1}{2}e^{-r^{2}}\right) \right\vert _{0}^{\infty }=\frac{\pi }{2}\left( \lim_{c\rightarrow \infty }-\frac{1}{2} e^{-c^{2}}+\frac{1}{2}e^{-0^{2}}\right) \\ &=&\frac{\pi }{2}\left( 0+\frac{1}{2}\right) =\frac{\pi }{4},\tag{5} \end{eqnarray*}$ yields $ I=\frac{\sqrt{\pi }}{2}.\tag{6} $