 Solution

In the translated coordinates $u=x-1$, $v=y+1$, the zero level set is the pair of lines $u=0$ and $v=0$. The function is positive in the quadrants $uv>0$ and negative in those with $uv<0$, and tends to zero as $u^2+v^2\to \infty$. Thus its nonzero contours form closed curves around the four off-centre stationary points.

For an explicit description of the contour plot, use polar coordinates $u=r\cos \theta$, $v=r\sin \theta$:At fixed direction within a quadrant, the magnitude increases until $r=\sqrt{2}$ and then decreases towards zero. Also $|u|e^{-u^2/2}\le e^{-1/2}$, with equality exactly when $|u|=1$, and similarly for $v$. Consequently $|f|\le e^{-1}$; each positive level below $e^{-1}$ encloses a positive peak in each positive quadrant, with the analogous negative contours in the negative quadrants.

To classify the stationary points, the Hessian matrix isAt $u=v=0$ it is $\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$, with eigenvalues $1$ and $-1$, so $(1,-1)$ is a saddle point. At $u,v\in \{-1,1\}$ it is $-2fI$. Therefore $(2,0)$ and $(0,-2)$ are strict local maxima, and $(2,-2)$ and $(0,0)$ are strict local minima. These are also global maxima and minima by the bound just proved.