Symmetry in double symsum from Symbolic math toolbox

Hello,

I am trying to implement the computation of moments for triangles, according to the quad sum given at

http://math.stackexchange.com/questions/3990/computing-the-moments-of-a-triangle

Unfortunately, I experience an asymmetric with respect to (x^p,y^q) pairs, i.e. for p=0 and q=2 I get the right expansion, while for p=2 and q=0, it seems incorrect.

I have used 4 times symsum separately, for easier read of the code, although I am not 100% sure at the substitution, whether the right indexes are matched.

As the code is rather short, I am am pasting here:

syms p q k l i j x0 x1 x2 y0 y1 y2 nu nusubs subsys moments;

nusubs = symsum( nchoosek(k,i) * nchoosek(l,j) / (k-i+l-j+1) * (x0-x1)^i * (x1-x2)^(k-i) * (y0-y1)^j * (y1-y2)^(l-j) , j, 0, l);

nu = symsum(nusubs, i, 0, k);

subsys = symsum( (-1)^(k+l) * nchoosek(p,k) * nchoosek(q,l)/(k+l+2) * nu * x0^(p-k) * y0^(q-l) , l, 0, q );

moments = symsum(2 * subsys, k, 0, p);

p = 2; q = 0;
subs(moments)
pretty(ans)

As I have no further idea for debugging this problem, I would highly appreciate any suggestion for improvement.

Regards, Levente