<s>[INST] Double Integral Question on unit square Hints on solving following double integral will be appreciated.
$$\int_0^1 \int_0^1 \frac{x^2-y^2}{(x^2+y^2)^2} \text{d}y \ \text{d}x$$
 [/INST] \begin{align}
\frac{x^{2} - y^{2}}{(x^{2}+y^{2})^{2}}
&=
{x^{2} + y^{2} -2y^{2} \over \left(x^{2} + y^{2}\right)^{2}}
=
{1 \over \left(x^{2} + y^{2}\right)^{2}}\,\left\lbrack%
{\partial y \over \partial y}\,\left(x^{2} + y^{2}\right)
-
{\partial\left(x^{2} + y^{2}\right) \over \partial y}\,y
\right\rbrack
\\[3mm]&=
{\partial \over \partial y}\left(y \over x^{2} + y^{2}\right)\,,
\qquad\qquad
\left(x, y\right) \not= \left(0, 0\right)
\\[5mm]&
\end{align}
\begin{align}
\lim_{\epsilon \to 0^{+}}\int_\epsilon^1 \int_0^1 \frac{x^2-y^2}{(x^2+y^2)^2} \text{d}y \ \text{d}x
&=
\lim_{\epsilon \to 0^{+}}\int_{\epsilon}^{1}\left.{y \over x^{2} + y^{2}}\right\vert_{\,y\ =\ 0}^{\,y\ =\ 1}\
{\rm d}x
=
\lim_{\epsilon \to 0^{+}}\int_{\epsilon}^{1}{{\rm d}x \over x^{2} + 1}
\\[3mm]&=
\arctan\left(1\right)
=
{\large{\pi \over 4}}
\end{align}
It seems to be the Mathematica package trick !!!: It's is equivalent to exclude $\left\lbrace \left(0,y\right)\ \ni\ y \in \left(0, 1\right) \right\rbrace$ and takes a limit after integration. It's like this:
\begin{align}
\int_{\epsilon}^{1}
\left.{y \over x^{2} + y^{2}}\right\vert_{\,y\ =\ 0}^{\,y\ =\ 1}\ {\rm d}x
=
\int_{\epsilon}^{1}{1 \over x^{2} + 1}\,{\rm d}x
=
{\pi \over 4} - \arctan\left(\epsilon\right) \to {\pi \over 4}
\end{align}
If you "exclude" both $x = 0$ and $y = 0$ you get:
$$
{\pi \over 4} - \arctan\left(\epsilon\right)-
\arctan\left(1 \over \epsilon\right) + \arctan\left(1\right) \to 0
$$
which is quite obvious ( the integral changes sign when $x \leftrightarrow y$ ).
The problem is related to the singular behavior of the 2D-Green function  of the Laplacian operator:
$\left.\nabla^{2}\ln\left(\rho\right)\right\vert_{\rho\ \not=\ 0} = 0$ but its integral around $\vec{\rho} = \vec{0}$ is $\not= 0$.
$\rho \equiv x\,\hat{x} + y\,\hat{y}$. It's common in 2D-Electrostatic.
For example 
\begin{align}
\int_0^1 \int_0^1 \frac{x^2-y^2}{(x^2+y^2)^2} \text{d}y \ \text{d}x
&=
\int_{S}\nabla\times\left({y \over x^{2} + y^{2}}\,\hat{x}\right)\cdot\hat{z}\,
{\rm d}x\,{\rm d}y
=
\oint\left({y \over x^{2} + y^{2}}\,\hat{x}\right)\cdot{\rm d}\vec{\rho}
\\[3mm]&=
\int_{1}^{0}{1 \over x^{2} + 1}\,\left(-{\rm d}x\right)
=
{\large{\pi \over 4}}
\end{align}
</s>