Double integral comes out different after changing $dx$ with $dy$.

I'm afraid your teacher is right. If we integrate over $x$ first,$$\int_0^1[(x+y)^{-2}-2y(x+y)^{-3}]dx=[y(x+y)^{-2}-(x+y)^{-1}]_0^1=-(1+y)^{-2}$$implies the result is $[(1+y)^{-1}]_0^1=-\frac12$. But if we integrate over $y$ first,$$\int_0^1[2x(x+y)^{-3}-(x+y)^{-2}]dy=[(x+y)^{-1}-x(x+y)^{-2}]_0^1=(x+1)^{-2}$$implies the result is $[-(1+x)^{-1}]_0^1=\frac12$.

Let's consider the geometry of this. The surface $z=\frac{x-y}{(x+y)^3}$ has parts with positive and negative $z$ for $(x,\,y)\in[0,\,1]^2$, and your hope was to define a net (above minus below) volume within this region for $x,\,y$. But the result is undefined, because the part-volumes above and below $z=0$ are each infinite, and $\infty-\infty$ is an indeterminate form. Simpler examples of this phenomenon (by which I mean they don't involve double integrals, but Wikipedia offers something similar that does) include $\int_{-a}^b\tfrac1xdx$ with $a,\,b>0$. As @CharlesHudgins notes, infinite series do this too.

If $\displaystyle \iint\limits_{[0,2]^2} \left| \frac{x-y}{(x+y)^3}\right| \, d(x,y) = +\infty$ then the two iterated integrals may have different (finite) values.

When you say, "this integral represents volume under a curve and the area over which we are integrating is same in both cases," the integral you are speaking of is a double integral. The two integrals you are computing are iterated integrals. When the function is absolutely integrable, then the double integral is equal either one of the iterated integrals, and of course, both iterated integrals have the same value. This is Fubini's theorem.

I haven't done the calculations here, but if what you say is true, it must be the case that the double integral is not absolutely integrable. In fact, it must not exist. What is the integral over the region $0\leq x \leq1, 0\leq y\leq x$?