Let $\gamma_1,\gamma_2:\mathbb{R}\rightarrow\mathbb{R}^2$ be defined by $\gamma_1(t) = (1+t,1+t)$ and $\gamma_2(t) = (1+t,1-t)$, and let $g_1 = f\circ\gamma_1$, $g_2 = f\circ\gamma_2$. Notice that $\gamma_1(0) = \gamma_2(0) = (1,-1)$, and
\begin{align} g_1(t) &= (1+t)+(1+t)-(1+t)(1+t) = 1-t^2 \\
g_2(t) &= (1+t)+(1-t)-(1+t)(1-t) = 1+t^2.
\end{align}
We thus see that $g_1$ attains a strict maximum at $t=0$, while $g_2$ attains a strict minimum at $t=0$. This means that $f$, restricted to the curve $\gamma_1$, attains a strict maximum at $(1,-1)$, while $f$ restricted to the curve $\gamma_2$ attains a strict minimum instead. Thus, $(1,-1)$ is a saddle point for $f$, as it is a critical point which cannot be a local maximum nor a local minimum.
How did I come up with the $\gamma_1,\gamma_2$? Basically, these are curves passing through the critical point $(1,1)$ at $t=0$ whose derivatives at $t=0$, namely $(1,1)$ and $(1,-1)$, are the eigenvectors of the Hessian $H(1,1) = \begin{pmatrix} 0 & -1 \\ -1 & 0 \end{pmatrix}$, of eigenvalues $-1$ and $1$, respectively. Why should these curves work? By Taylor's formula, we have
\begin{align} f(x_0+x)-f(x_0) &\approx (\nabla f(x_0))^Tx + \frac{1}{2}x^TH(x_0)x \\
&=\frac{1}{2}x^TH(x_0)x\quad\text{if }\nabla f(x_0) = 0\text{, i.e. }x_0\text{ is a critical point.}
\end{align}
Thus if $x = t\begin{pmatrix} 1 \\ 1 \end{pmatrix}$, then $x$ is an eigenvector of $H$ of eigenvalues $-1$, and hence
$$ \frac{1}{2}x^THx = \frac{1}{2}x^T(-x) = -\frac{1}{2}\left(t\begin{pmatrix} 1 \\ 1 \end{pmatrix}\right)^T\left(t\begin{pmatrix} 1 \\ 1 \end{pmatrix}\right) = -t^2$$
while similarly $x = t\begin{pmatrix} 1 \\ -1\end{pmatrix}\implies \frac{1}{2}x^THx = t^2$. Since $\gamma_1(t) = (1,1)+t(1,1)$ (and similarly for $\gamma_2$), by our arguments we have
$$ f(\gamma_1(t))-f(1,1) \approx -t^2 <0 $$
and
$$ f(\gamma_2(t))-f(1,1) \approx t^2 > 0. $$
To make the arguments fully rigorous, one has to replace the "$\approx$" in Taylor's formula by the actual error, which produces a term that vanishes faster than quadratically and hence will not affect our estimates. In this case the formula is actually exact, since we are just dealing with a second-order polynomial.