This one requires a trick. Let $I$ be the value of the integral. Then
$\begin{align*} I^2&=\left(\int_0^\infty\frac4{\pi b^2}x^2e^{-x^2/b^2}dx\right)\left(\int_0^\infty\frac4{\pi b^2}y^2e^{-y^2/b^2}dy\right)\\\\ &=\frac{16}{\pi^2b^4}\int_0^\infty\int_0^\infty x^2y^2e^{-(x^2+y^2)/b^2}dydx\;. \end{align*}$
Now convert to polar coordinates: $x=r\cos\theta$, $y=r\sin\theta$, etc. You’re integrating over the first quadrant, so you want your double integral in polar coordinates to have $0\le\theta\le\frac{\pi}2$ and $0\le r<\infty$. When you’ve completed the integration, you’ll have $I^2$, from which you can easily get $I$.
Added: Ignoring the various constants, you have essentially something like $\int_0^{\pi/2}\int_0^\infty r^5\cos^2\theta\sin^2\theta e^{-r^2}drd\theta=\int_0^{\pi/2}\cos^2\theta\sin^2\theta\int_0^\infty r^5e^{-r^2}drd\theta\;.$
The inner integral (with respect to $r$) can be done by repeated integration by parts; for the first one let $u=r^4$, $dv=re^{-r^2}dr$, so that $du=4r^3dr$ and $v=-\frac12e^{-r^2}$. That will leave you with something of the form $\int_0^\infty r^3e^{-r^2}dr$ to deal with. Repeat the process, and you’ll have something of the form $\int_0^\infty re^{-r^2}dr$, which you can integrate outright.
At that point you’ll be integrating some multiple of $\cos^2\theta\sin^2\theta$. One way is to use the double angle formula for the sine to rewrite this as $\frac12\sin^22\theta$, then use the half-angle formula to rewrite $\sin^22\theta$ as $\frac12(1-\cos 4\theta)$, which you can integrate.