In[]:=
MultivarTaylor3D[f_,vars_,point_,order_]:=Module{L=order,X1=vars[[1]],X2=vars[[2]],X3=vars[[3]],P1=point[[1]],P2=point[[2]],P3=point[[3]]},TermList=Array(D[f,{X1,(#1-1)},{X2,(#2-1)},{X3,(#3-1)}]/.{X1P1,X2P2,X3P3})&,{L,L,L}//FullSimplify;Return[TermList];
1
(#1-1)!(#2-1)!(#3-1)!
(#1-1)
(X1-P1)
(#2-1)
(X2-P2)
(#3-1)
(X3-P3)
In[]:=
TestFunction=;
3
(x+1)
3
(y+1)
3
(z+1)
In[]:=
TestArray=MultivarTaylor3D[TestFunction,{x,y,z},{0,0,0},3]//FullSimplify
Out[]=
{{{1,3z,3},{3y,9yz,9y},{3,9z,9}},{{3x,9xz,9x},{9xy,27xyz,27xy},{9x,27xz,27x}},{{3,9z,9},{9y,27yz,27y},{9,27z,27}}}
2
z
2
z
2
y
2
y
2
y
2
z
2
z
2
z
2
y
2
y
2
y
2
z
2
x
2
x
2
x
2
z
2
x
2
x
2
x
2
z
2
x
2
y
2
x
2
y
2
x
2
y
2
z