I honestly don't understand why your subs call doesn't work. I can't make it work either, in spite of a number of different approaches that have worked for me previously. I consider this a bug, so you may want to report it to TMW as such.
All is not lost, however. Here's a kludgy but effective work-around. Note that it still only recognizes (dx*dy) and not such as (dx*dy)^2:
syms dx dy x y
PO = expand( (x + dx) * (y + dy) + (dx * dy) ) % Prototype expression
[POc, T] = coeffs(PO, [dx dy])
Tfind = find(T == 'dx*dy')
POc(Tfind) = 0
PO = POc * T.'
This could probably be made more efficient. I wrote it as I did to make the logic behind it transparent.
