now f(0,b)=0 and f(a,0)=0 thus we have

f(a,b)=k*ab*g(a,b)

where k is a constant and g(a,b) is a polynomial in a,b

now if a=-b then f(a,b)=0 as well thus we can simplify further to

f(a,b)=k*ab*(a+b)*g(a,b)

again with k a constant and g(a,b) a polynomial in a,b

now by expanding f(a,b) and dividing out ab(a+b) we can see that

k=7 and g(a,b)=a^4+2a^3b+3a^2b^2+2ab^3+b^4

now this is where a bit of intuition comes into play.  The symmetry in the coefficients makes me think of squared polynomial thus we have

g(a,b)=[h(a,b)]^2

now we can deduce that h(a,b) is a polynomial in a,b of degree two with 3 terms two of which are a^2 and b^2 and since the remaining terms of g(a,b) have a common factor of ab then the third term of h(a,b) must be ab

thus we can check if

g(a,b)=(a^2+ab+b^2)^2

which by expanding is verified

thus we are left with

(a+b)^7-a^7-b^7=7ab(a+b)(a^2+ab+b^2)^2