This old code
Remove[bilgicTaylor];
bilgicTaylor[f_, varL_List?VectorQ, zeroL_List?VectorQ,
degree_Integer?NonNegative] :=
Module[{X, tPD},
(* total preDifferential *)
X /: X[o_List]^n_Integer := Product[X[o], {i0, n}];
X /: X[o_List] X[p_List] := X[o + p];
tPD = Plus @@ (zeroL*
Thread[X[
Sort[Permutations[Join[{1}, Table[0, {Length[zeroL] - 1}]]],
Flatten[Position[#, 1] &]]]]);
(* result *)
f[Sequence @@ varL] + Sum[Expand[tPD^d]/d!, {d, 1, degree}] //.
X[oo_List] :> Derivative[Sequence @@ oo][f][Sequence @@ varL]
] /; Length[varL] == Length[zeroL]
should do it
In[18]:= Clear[k]
k[u_, v_] := (u + v)^2
In[20]:= bilgicTaylor[k, {x, y}, {x0, y0}, 2]
Out[20]= 2 x0 (x + y) + (x + y)^2 + 2 (x + y) y0 + 1/2 (2 x0^2 + 4 x0 y0 + 2 y0^2)
In[23]:= %20 // Simplify
Out[23]= (x + x0 + y + y0)^2
check it simply by undefining the k
Remove[k]
bilgicTaylor[k, {x, y}, {x0, y0}, 2]
