This is a naive string manipulation, but it works OK for positive integer powers.

    newTerm[a_,b_]:=ToExpression[(a//ToString)<>"¢"<>(b//ToString)]

    collect[expr_,{vars__}]:=StringReplace[ToString[Collect[expr/.Power->newTerm,Flatten[{{vars}/.Power->newTerm}]]],"¢"->"^"]//ToExpression

    collect[expr_,vars_]:=StringReplace[ToString[Collect[expr/.Power->newTerm,Flatten[{{vars}/.Power->newTerm}]]],"¢"->"^"]//ToExpression

For example:

In[179]:= collect[a x+b x+c x^2+d x^2+e y + f y+g y^2+h y^2,{x^2,y^2,x}]
Out[179]= (a+b) x+(c+d) x^2+e y+f y+(g+h) y^2

or, with x^2 y collected:

In[180]:= collect[a x+b x+c x^2 y+d x^2 y+e y + f y+g y^2+h y^2,{x^2 y,y^2,x}]
Out[180]= (a+b) x+e y+f y+(c+d) x^2 y+(g+h) y^2

or, with only one term to collect, the braces are not needed:

In[211]:= collect[a x^2+b x+c x^2+d x,x^2]
Out[211]= b x+d x+(a+c) x^2

With more variable you have to decide how to deal with mixed terms x1*x2

A similar problem was discussed in the thread "How to Get a Taylor Series for Multiple Variables?".

Thank you, enough straightforward. Can the code be extended for expression with several variables (i.e., x1^2, x2^2, ..., xn^2)?