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