How to prove that integral $$\iint\limits_{x\ge1,~y\ge1}\frac{x^2-y^2}{(x^2+y^2)^2}\,dx\,dy$$ is divergent?
Let $x=1+r\cos\varphi$, $y=1+r\sin\varphi$, then $$\iint\limits_{x\ge1,~y\ge1}\dfrac{x^2-y^2}{(x^2+y^2)^2}\,dxdy=\iint\limits_{\substack{0\le\varphi\le\frac{\pi}2\\r\ge0}}\dfrac{r^3\cos2\varphi+2r^2(\cos\varphi-\sin\varphi)}{(r^2+2r(\cos\varphi+\sin\varphi)+2)^2}\,drd\varphi.$$ Because $$\dfrac{|r^3\cos2\varphi+2r^2(\cos\varphi-\sin\varphi)|}{(r^2+2r(\cos\varphi+\sin\varphi)+2)^2}\ge\dfrac{|r^3\cos2\varphi-2r^2|} {(r^2+4r+2)^2}\ge\dfrac{r^3|\cos2\varphi|-2r^2}{(r^2+4r+2)^2},$$ we get that the integral $$\iint\limits_{x\ge1,~y\ge1}\frac{x^2-y^2}{(x^2+y^2)^2}\,dx\,dy$$ will diverge if the integral diverges \begin{multline*} \iint\limits_{\substack{0\le\varphi\le\frac{\pi}2\\r\ge0}}\dfrac{r^3|\cos2\varphi|-2r^2}{(r^2+4r+2)^2}\,drd\varphi=\\ =\iint\limits_{\substack{0\le\varphi\le\frac{\pi}2\\r\ge0}}\dfrac{r^3|\cos2\varphi|}{(r^2+4r+2)^2}\,drd\varphi-\iint\limits_{\substack{0\le\varphi\le\frac{\pi}2\\r\ge0}}\dfrac{2r^2}{(r^2+4r+2)^2}\,drd\varphi=\\ =\int\limits_0^{\pi/2}|\cos2\varphi|\,d\varphi\int\limits_0^{+\infty}\dfrac{r^3\,dr}{(r^2+4r+2)^2}-\int\limits_0^{\pi/2}d\varphi\int\limits_0^{+\infty}\dfrac{2r^2\,dr}{(r^2+4r+2)^2}=\\ =\int\limits_0^{+\infty}\dfrac{r^3\,dr}{(r^2+4r+2)^2}-\frac{\pi}2\int\limits_0^{+\infty}\dfrac{2r^2\,dr}{(r^2+4r+2)^2}. \end{multline*} What follows is obvious.
The integral does not exist (converge) in the sense of Lebesgue, since
$$\int_{[1,\infty)^2} \left|\frac{x^2-y^2}{(x^2+ y^2)^2} \right| = + \infty$$
The correct way to prove this when changing to polar coordinates is to consider the absolute value of the integrand. Since $\cos 2\varphi$ changes sign over the interval $[0,\pi/2]$, there is sufficient cancellation for the iterated improper integral to converge.
However, we can show that the iterated integrals converge directly without changing variables-- albeit to values with different signs depending on the order.
Note that
$$\frac {x^2-y^2}{(x^2+y^2)^2} = \frac{\partial}{\partial y} \left(\frac{y}{x^2 + y^2}\right) = -\frac{\partial}{\partial x} \left(\frac{x}{x^2 + y^2}\right), $$
and so
$$\begin{align}\int_1^\infty \left(\int_1^\infty \frac {x^2-y^2}{(x^2+y^2)^2}\, dx\right)\, dy &= \int_1^\infty \left(\int_1^\infty -\frac{\partial}{\partial x}\left(\frac {x}{x^2 + y^2}\right)\, dx\right)\, dy \\ &= \int_1^\infty \left.\frac{-x}{x^2+y^2}\right|_{x = 1}^{x = \infty} \, dy\\ &= \int_1^\infty \frac{dy}{1+y^2} \\ &= \frac{\pi}{4}\end{align}$$
By interchanging $x$ and $y$ it is easy to see that the iterated integral in reverse order takes the value $- \pi/4$.
Absolute divergence
We have
$$\int_{[1,\infty)^2}\left|\frac{x^2-y^2}{(x^2+ y^2)^2} \right| \geqslant \int_0^{\pi/2} \int_{\sqrt{2}}^R \frac{|\cos 2 \varphi|}{r} \, dr \, d \varphi = \int_0^{\pi/2}|\cos 2 \varphi| \, d \varphi \int_{\sqrt{2}}^R \frac{dr}{r}$$
and the RHS diverges as $R \to \infty$. Note that the iteration of integrals on the RHS is independent of order as the integrand is continuous and the region of integration is bounded.