I would like to compute the Taylor expansion of monomials whose powers are nonnegative integers variables.

First, here is a function which encodes the multivariate Taylor expansion at the points pts, at the order ord and with respect to the variables var.

mTaylor[expr_,var_List,pts_List,ord_Integer] /; Length[var] == Length[pts] := \
    Normal[(expr /. Thread[var -> s (var-pts)+pts]) + O[s]^ord] /. {s->1}

Assume I want to compute the multivariate Taylor expansion of the expression:

expr = x^a1 + x^a2 y^3

at the order $2$, wrt $(x,y)$ at $(0,0)$ and where $a1,a2$ are nonnegative integers.

The expected results would be

res = x^a1 (*=mTaylor[expr,{x,y},{0,0},3]*)

However I can't get my Taylor function correct to yield this result.

I tried using global assumptions :

$Assumptions = Element[a1,NonNegativeIntegers] && Element[a2,NonNegativeIntegers]

I tried also to write a routine to remove the high order terms but again it doesn't seem it handles the assumptions on the power:

EliminTe[expr_,var_List,ordre_Integer]:=Module[{},\
      FromCoefficientRules[Select[CoefficientRules[expr, var], Refine[Total@#[[1]] <= ordre    ,Assumptions->Element[a2,NonNegativeIntegers]] &], var]] 

Here is my code

mTaylor[expr_,var_List,pts_List,ord_Integer] /; Length[var] == Length[pts] := \
Normal[(expr /. Thread[var -> s (var-pts)+pts]) + O[s]^ord] /. {s->1}

$Assumptions = Element[a1,NonNegativeIntegers] && Element[a2,NonNegativeIntegers]

expr = x^a1 + x^a2 y^3 

res = Refine[mTaylor[expr,{x,y},{0,0},3]]