I will assume that $f$ is "nice". This isn't comprehensive; but will allow you to find, in most cases:
1) The local extrema of $f$ in the interior of its domain.
2) The global extrema of $f$ if $f$ is defined on a "nice" region $R$ with boundary $B$.
To find local extrema in the interior of the domain of $f$:
First, find the points where the gradient of $f$ is the zero vector or does not exist. These are the candidates (and only candidates) for local extreme values in the interior of the domain of $f$.
Given one these points, the so-called "second partials test" might be used to determine if the point gives a local extreme value. I say "might" since, in some cases, this test gives no useful info. If this turns out to be the case, then you need to examine what happens near the critical point in question by looking at $f$ itself (this may be hard).
To find global extrema where $f$ is defined on a "nice" region $R$ with boundary $B$:
You will need to examine the function on the boundary. Here, you will need to find a function of one variable, $g$, that gives the same value as $f$ at any point on the boundary. Usually, you find a parameterization of $B$ , and then substitute into $f$ to obtain $g$. Then find the global extrema of $g$.
You then compare these values with the values of $f$ that are candidates for local extrema in the interior.
Evaluate $f$ at each point in the interior where the gradient is the zero vector or does not exist. Compare these values with the global extreme values found on the boundary and select the largest and smallest of these values.
Candidate points ${\bf x}=(x,y)$ for local extrema:
A function of two variables $f$ can assume an extreme value only at points where
$
\nabla f({\bf x})=\bf 0
$ or does not exist. Here, $\nabla f$ is the gradient of $f$:
$$
\nabla f({\bf x})=f_x({\bf x})\,{\bf i}+f_y({\bf x})\,{\bf j}.
$$
We call any point where $
\nabla f({\bf x})=\bf 0
$ or does not exist a critical point of $f$.
Second partials test:
Suppose $f$ has continuous second partial derivatives.
For a critical point $ {\bf x}=(x,y)$ of $f$ in the interior of the domain of $f$,
set $$
D=f_{xx}({\bf x})f_{yy}({\bf x}) - f_{xy}^2({\bf x}).
$$
Then
$$
\eqalignno{
D\gt0 \hbox{ and } f_{xx}({\bf x})\gt0 \quad & \Longrightarrow \quad f({\bf x})\ \hbox{is a local min. value.} \cr
D\gt0 \hbox{ and }f_{xx}({\bf x})\lt 0 \quad & \Longrightarrow\quad f({\bf x})\ \hbox{is a local max. value.}\cr
D\lt 0 \quad & \Longrightarrow\quad f({\bf x})\ \hbox{is a saddle point.}\cr
}
$$
Finding Global Extrema
1) Find the critical points of $f$ contained in the interior of its domain. Evaluate $f$ at each of these points.
2) Find the global maximum and minimum value of $f$ on the boundary of the domain of $f$.
If the domain of $f$ is all of $\Bbb R^2$, skip this step.
3) Compare the values found in 1) and 2) and select the largest and smallest of these to obtain, respectively, the global maximun and minimum value of $f$.
Putting $f_x=0$ and $f_y=0$ you obtained the two equations $yx^2=1$ and $xy^2=1$. The solution of these two equations is not the line $x=y$, as you seem to think, but the single point $p:=(1,1)$. This is the only critical point in the domain ${\mathbb R}^2\setminus\{(x,y)\ |\ xy=0\}$ of $f$. One computes
$$f_{xx}(p)=2, \quad f_{xy}(p)=1,\quad f_{yy}(p)=2\ ,$$
whence $f_{xx}(p)>0$ and $D(p):=(f_{xx}f_{yy}-f_{xy}^2)_p=3>0$. It follows that $f$ takes a local minimum at $p$.
When you restrict $f$ to the open first quadrant then $f$ in fact takes its global minimum value $3$ at $p$. This can be seen as follows:
$$xy+{1\over x}+{1\over y}=\sqrt{y}\Bigl(x\sqrt{y}+{1\over x\sqrt{y}}\Bigr)+{1\over y}\geq 2\sqrt{y}+{1\over y}\qquad(x>0)\ ,$$
and here the right side has minimal value $3$ for $y=1$.
------- (addendum for extra information) ---------
The last step requires
$2\sqrt{y}+{1\over y} \ge 3$ for $y >0$. To see this,
write $y=(1+z)^2$ with $z > -1$. Then we want to establish
$$
2{(1+z)} + \frac{1}{(1+z)^2} \ge 3 \\
\leftrightarrow
2{(1+z)^3} + 1 - 3 {(1+z)^2} \ge 0\\
\leftrightarrow
z^2(2z +3) \ge 0
$$
which is true since $z^2 \ge 0$ and since for $z > -1$, $2z +3 > 1$. Also, the inequality is tight since equality occurs at $z=0$, i.e. $y=1$.
Best Answer
Hint: Look at the corresponding closed region $R$, so that the extrema will either be critical points in int$R$ or they will lie on the boundary of $R$. For the saddle points, use the second derivative test.
Second hint:
Step 1:
Solve the following equations and use the points in the second derivative test to decide if they are local extrema or saddle points. Suggestion: let $z=y/x$
$f_{x}=4x^{3}-4y=0\\f_{y}=4y^{3}-4x=0$
Step 2:
Now consider each piece of the boundary, making the appropriate substitution and extremize the function (of a single variable) you get in each case. For example, if we look at the top side of the square, where $y=2$ and $-2<x<2$, we get
$f(x,2)=x^{4}-8x+16$, from which $f'(x,2)=4x^{3}-8=0$ gives $x=2^{1/3}$ as a critical point. Now determine what type of critical point it is, substitute into $f(x,y)$ and compare it with what you got in step 1. Then check $x=2, x=-2$ in $f(x,2)$ and write down the result.
repeat for the other sides of the square.
All told, you get a table of values from which you can read off the maxima and minima.