THIS IS A PARTIAL ANSWER
To not have to deal with the difficult notation of limits (having to write the whole limit everytime) in order to prove that the function f(x,y) is differentiable at a fixed point (x,y), it is sufficient to write f(x+h,y+k)-f(x,y) as grad(f) . (h,k) + O(h,k) meaning f(x+h,y+k)-f(x,y) = grad(f) . (h,k) + O(h,k)
Where . is a scalar product and O(h,k) is something that tends to zero faster than (h,k). This approach is equivalent (and way easier to understand) to the limit notation.
You just that take your function f(x,y) = xy(x+y). You do as I did in the figure. The main thing is to expand everything and after, put together in ont  the terms with variable h^1 only (not k, and not h^2...) and the ones with k^1 only. These will be our scalar product grad(f) . (h,k) since it must be like ...h+...k. We put the rest of the terms together in an other . It will be the error of approximation as we call it : the O(h,k). We have to prove that this O(h,k) tends to zero faster than (h,k). After calculations, I find that it does work (you see the limit = 0 because in the nominator the sum of the powers of h and k are equal or higher than 2 in every term so it goes faster to zero than the h and k in the denominator which only have power 1 i.e. : lim (x->0) x^2/x = lim (x->0) x =0).
So for a fixed point (x,y) we succeeded in writing f(x+h,y+k)-f(x,y) = grad(f) . (h,k) + O(h,k) where O(h,k) tends faster to 0 than (h,k).
I know it is not easy compared to function with one variable but going with the limits is even more difficult in notation.
I hope it helps.
I do not want to complete answering all your questions until you understand this concept. So if you understand this just post back that you understand and that you want us to continue answering the reste of the questions. If not, tell us you do not understand.
In any case, I hope that me or one of my collegues will help you.