Wolfgang's solution is much better than mine (when it is applicable).
There is one thing to watch out for, though. Maxima simplifies (e.g.)
x+(a+b) to x+a+b, so facsum can't collect the terms in x^0*y^0.
But there is a simple workaround:
subst([x=1,y=1,foo=1], args( facsum( foo * expr, x, y, foo)))$
-s