Solution (source code)

= Solution
{c}

Write $A=a^2>0$ and $E=e^{-(x^2+y^2)/A}>0$. The <partial derivatives> factor as
$$
f_x=2xE\left(1-\frac{x^2-2y^2}{A}\right),\qquad
f_y=-2yE\left(2+\frac{x^2-2y^2}{A}\right).
$$
A <critical point> has both <partial derivatives> zero. If $x,y\ne0$, this would require $x^2-2y^2=A$ and $x^2-2y^2=-2A$, which is impossible. On the axes the remaining conditions give
$$
\boxed{(0,0),\quad(\pm|a|,0),\quad(0,\pm|a|).}
$$
At either <critical point> with $x\ne0$, the <Hessian matrix> is
$$
H_f=\begin{pmatrix}-4/e&0\\0&-6/e\end{pmatrix}.
$$
It is <negative definite>, so the <second-derivative test> shows that \b[both $(\pm|a|,0)$ are strict <local maxima>], of value $a^2/e$. They also attain the <global maximum>: $f\le x^2e^{-x^2/A}\le A/e$. For completeness, the <Hessian matrix> at $(0,\pm|a|)$ is $\operatorname{diag}(6/e,8/e)$, giving strict <local minima>, while at the origin it is $\operatorname{diag}(2,-4)$, giving a <saddle point>.