Solution (source code)

= Solution
{c}

For an unconstrained differentiable <function> of two real variables, an interior <local minimum> must solve the simultaneous <stationary point> equations
$$
f_x(x,y)=0,\qquad f_y(x,y)=0,
\quad\text{or}\quad \nabla f(x,y)=0.
$$
Apply <Newton's method in optimization> to this <gradient> equation. The <Jacobian matrix> of the <gradient> is the <Hessian matrix>
$$
H(x,y)=\begin{pmatrix}f_{xx}&f_{xy}\\f_{yx}&f_{yy}\end{pmatrix}.
$$
At $p_n=(x_n,y_n)^T$, solve the <linear system> $H(p_n)s_n=-\nabla f(p_n)$, then set $p_{n+1}=p_n+s_n$. When the <Hessian matrix> is invertible, this is
$$
\boxed{p_{n+1}=p_n-H(p_n)^{-1}\nabla f(p_n).}
$$
This follows by setting the first order <Taylor expansion> $\nabla f(p_n+s_n)\simeq\nabla f(p_n)+H(p_n)s_n$ equal to zero. Solving the <linear system> avoids explicitly computing a <matrix inverse>. A <stationary point> need not be a <local minimum>: for example, a <positive-definite matrix> $H(p)$ suffices to classify a <stationary point> as a strict <local minimum>, whereas an indefinite <Hessian matrix> identifies a <saddle point of a scalar function>.