I would like to either confirm my solution is correct, or find the error in it. I used the following MATLAB code to try to check my answer, but the solution it gave differs from mine.
f = exp(-(x^2+y^2))cos(xy);
taylor(f,[x,y],[0,0])
ans =
x^4/2 + (x^2*y^2)/2 - x^2 + y^4/2 - y^2 + 1
My solution:
The second degree Taylor polynomial at the point $(a,b)$ is given by $$p_2(x,y)=f(a,b)+Df(a,b)\left[\begin{array}{c} x-a\\y-b \end{array}\right]+\frac{1}{2}[x-a\mbox{ }y-b]Hf(a,b)\left[\begin{array}{c} x-a\\y-b \end{array}\right]$$ where $Df(a,b)$ is the Jacobian matrix and $Hf(a,b)$ is the Hessian matrix (i.e. first and second order partial derivatives, respectively). Thus, we begin by computing the partial derivatives: \begin{equation*} \begin{split} \frac{\partial f}{\partial x}(x,y)=-e^{-x^2-y^2}(y\mbox{ sin}(xy)+2x\mbox{ cos}(xy))\\ \frac{\partial f}{\partial y}(x,y)=-e^{-x^2-y^2}(x\mbox{ sin}(xy)+2y\mbox{ cos}(xy))\\ \frac{\partial^2 f}{\partial x^2}(x,y)=e^{-x^2-y^2}((4x^2-y^2-2)\mbox{cos}(xy)+4xy\mbox{ sin}(xy))\\ \frac{\partial^2 f}{\partial y^2}(x,y)=e^{-x^2-y^2}(4xy\mbox{ sin}(xy)-(x^2-4y^2+2)\mbox{cos}(xy))\\ \frac{\partial^2 f}{\partial x\partial y}(x,y)=\frac{\partial^2 f}{\partial y\partial x}(x,y)=e^{-x^2-y^2}((2x^2+2y^2-1)\mbox{sin}(xy)+3xy\mbox{ cos}(xy)) \end{split} \end{equation*} Then, at the point $(x_0,y_0)=(0,0)$, \begin{equation*} \begin{split} f(0,0)=1\\ \frac{\partial f}{\partial x}(0,0)=\frac{\partial f}{\partial y}(0,0)=0\\ \frac{\partial^2 f}{\partial x^2}(0,0)=\frac{\partial^2 f}{\partial y^2}(0,0)=-2\\ \frac{\partial^2 f}{\partial x\partial y}(0,0)=\frac{\partial^2 f}{\partial y\partial x}(0,0)=0 \end{split} \end{equation*} Therefore, the second degree Taylor polynomial at the point $(0,0)$ is \begin{equation*} \begin{split} p_2(x,y)=f(0,0)+Df(0,0)\left[\begin{array}{c} x-0\\y-0 \end{array}\right]+\frac{1}{2}[x-0\mbox{ }y-0]Hf(0,0)\left[ \begin{array}{c} x-0\\y-0 \end{array}\right]\\ =1+[0\mbox{ }0]\left[ \begin{array}{c} x\\y \end{array}\right]+\frac{1}{2}[x\mbox{ }y]\begin{bmatrix} -2&0\\0&-2 \end{bmatrix}\left[ \begin{array}{c} x\\y \end{array}\right]=1-x^2-y^2 \end{split} \end{equation*}
I leave your answer's critique to someone with less tired eyes. I will tell you this much. There is an easier way to calculate: just substitute into known Taylor series for the exponential and cosine, $$f(x,y) = e^{-x^2-y^2}\cos(xy) = (1-x^2-y^2+\cdots)(1-\frac{1}{2}(xy)^2+ \cdots )$$ But, $(xy)^2$ is order $4$ so the answer is merely, $$f(x,y) = e^{-x^2-y^2}\cos(xy) = 1-x^2-y^2+\cdots $$ The above is the Taylor expansion of $f$ at $(0,0)$ to order $2$.