Here is the function:

$$f(x,y)=e^{xy}sin(x^2y), \ \ M=\left[\begin{matrix}1 \\ \pi \end{matrix}\right], n=4, k=3$$

$$f(x,y)= {\Sigma_{n=0}^{\infty}} \left[ \left(\begin{matrix} n \\ k \end{matrix}\right)\left. \frac{\partial f}{ \partial x^{n-k}\partial y^k }\right\vert_{M} (x-x_0)^{n-k}(y-y_0)^k \right]$$

My solution:

$$f_x=ye^{xy}sin(x^2y) + 2xe^{xy}cos(x^2y)$$

$$f_{xy}=e^{xy}sin(x^2y)+xye^{xy}sin(x^2y)+ x^2ye^{xy}cos(x^2y)+2x^2e^{xy}cos(x^2y)-2x^3e^{xy}sin(x^2y)$$

$$f_{xyy}=xe^{xy}sin(x^2y)+x^2e^{xy}cos(x^2y)+xe^{xy}sin(x^2y)+x^2ye^{xy}sin(x^2y)+x^3ye^{xy}cos(x^2y)+x^2e^{xy}cos(x^2y)+x^3ye^{xy}cos(x^2y)-x^4ye^{xy}sin(x^2y)+2x^3e^{xy}cos(x^2y)-2x^4e^{xy}sin(x^2y)-2x^4e^{xy}sin(x^2y)-2x^5e^{xy}cos(x^2y)$$

$$f_{xyyy}=x^2e^{xy}sin(x^2y)+x^3e^{xy}cos(x^2y)+x^3e^{xy}cos(x^2y)-x^4e^{xy}sin(x^2y)+x^2e^{xy}sin(x^2y)+x^3e^{xy}cos(x^2y)+x^2e^{xy}sin(x^2y)+x^3ye^{xy}sin(x^2y)+x^4ye^{xy}cos(x^2y)+x^3e^{xy}cos(x^2y)+x^4ye^{xy}cos(x^2y)-x^5ye^{xy}sin(x^2y)+x^3e^{xy}cos(x^2y)-x^4e^{xy}sin(x^2y)+x^3e^{xy}cos(x^2y)+x^4ye^{xy}cos(x^2y)-x^5ye^{xy}sin(x^2y)-x^4e^{xy}sin(x^2y)-x^5ye^{xy}sin(x^2y)-x^6ye^{xy}cos(x^2y)+2x^4e^{xy}cos(x^2y)-2x^5e^{xy}sin(x^2y)-2x^4e^{xy}sin(x^2y)-2x^6e^{xy}cos(x^2y)-2x^5e^{xy}sin(x^2y)-2x^6e^{xy}cos(x^2y)-2x^6e^{xy}cos(x^2y)+2x^7e^{xy}sin(x^2y)$$

So, I know that it is possible to write it in a shorter way, but I don't know how. Because for me it looks so messy the way I did it. So then: $$f(1,\pi)= \frac{4!}{3!}\cdot( e^{\pi}+ e^{\pi}+e^{\pi}+\pi e^{\pi}+e^{\pi}+\pi e^{\pi}+e^{\pi}+e^{\pi}+\pi e^{\pi}-\pi e^{\pi}+2e^{\pi}-2e^{\pi}-2e^{\pi})= 6 \cdot (4\pi +\pi e^{\pi})$$

So, the Taylor expansions term that corresponds to given $n$ and $k$ is: $$6\cdot(4\pi + \pi e^{\pi})(x-1)(x-\pi)^3$$

Any feedback would be helpful.