collect multiple variables and simplify their coefficients

edit

maxima allows to collect multiple terms, so it is possible to do so in Sage using:

sage: f = 4*x*y + x*z + 20*y^2 + 21*y*z + 4*z^2 + x^2*y^2*z^2
sage: fm = f.maxima_methods()
sage: fm.collectterms(x,y)
x^2*y^2*z^2 + 4*x*y + x*z + 20*y^2 + 21*y*z + 4*z^2
sage: fm.collectterms(x)
x^2*y^2*z^2 + (4*y + z)*x + 20*y^2 + 21*y*z + 4*z^2

from there you could use the operands() method to simplify each term separately. Is this enough?

Can you assume that your expression is polynomial in x and y?

Here are two helpful functions:

sage: var('x,y,z')
(x, y, z)
sage: expr=x^2*y+2*x*y*z+54*z-3*x*z
sage: expr.collect(x)
x^2*y + (2*y*z - 3*z)*x + 54*z
sage: expr.coefficients(x)
[[54*z, 0], [2*y*z - 3*z, 1], [y, 2]]