Find the second degree Taylor polynomial of $f(x,y)=e^{-x^2-y^2}cos(xy)$ at $x_0=0$, $y_0=0$. 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$.
