Let $f:\mathbb{R}^2 \mapsto \mathbb{R}^2$. If $f$ is twice derivable, then for $h=\left(h_1,h_2\right)$, the Taylor polynomial in point $a$ is:
$$ f\left(a+h\right)= f\left(a\right) + h_1\frac{\partial f}{\partial x_1}\left(a\right) + h_2\frac{\partial f}{\partial x_2} + \frac{1}{2}\left(h_1^{2}\frac{\partial^2 f}{\partial x_1^{2}}\left(a\right) + 2h_1h_2\frac{\partial^2f}{\partial x_1x_2}\left(a\right) + h_2^2\frac{\partial^2 f}{\partial x_2^2}\left(a\right)\right) + \underset{h→ 0}{o}\left(||h||^2\right) $$
So taking this in mind, I tried to calculate the Taylor Polynomial of $f:\left(x,y\right)\mapsto \left(ye^x,\cos \left(xy\right)\right)$.
So I got this:
$$\frac{\partial f}{\partial x}\left(x,y\right)=\left(ye^x,-y\sin \left(xy\right)\right)$$
$$\frac{\partial f}{\partial y}\left(x,y\right)=\left(e^x,-x\sin \left(xy\right)\right)$$
$$\frac{\partial^2 f}{\partial x^2}\left(x,y\right)=\left(ye^x,-y^2\cos \left(xy\right)\right)$$
$$\frac{\partial^2 f}{\partial y^2}\left(x,y\right)=\left(1,-x^2\cos \left(xy\right)\right)$$
$$\frac{\partial^2}{\partial x \partial y}\left(x,y\right)=\left(e^x,-\sin \left(xy\right)-yx\cos \left(xy\right)\right) $$
$$\frac{\partial^2}{\partial y \partial x}\left(x,y\right)=\left(e^x,-\sin \left(xy\right)-yx\cos \left(xy\right)\right) $$(the same based on the continuity and Schwartz theorem)
So for the Taylor Polynomial at $0$, do we get this? ( MSE seems to have troubles rendering my Latex code, so I made a screenshot of it. The source code will be in a link below)
Source code: https://justpaste.it/11pgy
On the right hand side, the function and the derivatives are evaluated at the point $a$, which is $(0,0)$ in this case. Hence $$f(x,y) \sim f(0,0) + x\frac{\partial f}{\partial x}(0,0) + y\frac{\partial f}{\partial y}(0,0) + \frac{1}{2} \left(x^2 \frac{\partial^2 f}{\partial x^2}(0,0) + 2xy\frac{\partial^2 f}{\partial x \partial y}(0,0) + y^2 \frac{\partial^2 f}{\partial y^2}(0,0) \right)$$
The Taylor polynomial of order two is the right hand side of the above and does not include the $o()$ part. Note that the output of $f$ and of its partial derivatives are in $\mathbb{R}^2$, hence the result will be a vector with components polynomials of degree at most two in $x$ and $y$:
$$(0,1) + x(0,0) + y(1,0) + \frac{1}{2} \left(x^2 (0,0) + 2xy (1,0) + y^2(0,0) \right) = (y+xy, 1) $$
[Note that the first component of $\frac{\partial^2 f}{\partial y^2}(x,y)$ is not 1.]
You could also use the Taylor/Maclaurin expansions for $e^u = 1+u + \frac{1}{2}u^2 + \dotsb $ and $\cos{u} = 1 - \frac{1}{2} u^2 + \dotsb$ to see that $$ye^x = y+xy + \text{terms of total degree at least 3}$$ and $$\cos{xy} = 1 +\text{terms of total degree at least 3}.$$