I am trying to find maximum of
\begin{equation} f(x, y) = x^2 - xy + y - 4x \end{equation}
\begin{equation}\label{constraints} \text{s.t. } 0 \leq x \leq 2 \text{ and } 0 \leq y \leq 1 \end{equation}
I got an advice to first check for the interior solution (it has none there is one saddle point at [1,-2] and afterwards check each constraint separately, but I am not sure if what I did is correct and how to proceed with finding extreme points at the boundaries
\begin{equation} \begin{aligned} 0=x; 0=y \\ 0=x; 1=y \\ 2=x; 0=y \\ 2=x; 1=y \\ \end{aligned} \end{equation}
Now we can set up 4 Lagrangians to check for each boundary:
\subsubsection*{First Lagrangian:}
\begin{equation} \mathcal{L} = x^2 - xy + y - 4x - \lambda_1(-x) - \lambda_2 (-y) \end{equation}
with FOCs
\begin{equation} \begin{aligned} \frac{\partial \mathcal{L} }{\partial x} & = 2x - y -4 + \lambda_1 = 0\\ \frac{\partial \mathcal{L}}{ \partial y} & = -x + 1 + \lambda_2 = 0 \\ x&=0\\ y&=0 \end{aligned} \end{equation}
hence here $x=y=0$ and $\lambda_1 = 4$ and $\lambda_2=-1$
\subsubsection*{Second Lagrangian:}
\begin{equation} \mathcal{L} = x^2 - xy + y - 4x - \lambda_1(-x) - \lambda_2 (-y) \end{equation}
with FOCs
\begin{equation} \begin{aligned} \frac{\partial \mathcal{L} }{\partial x} & = 2x - y -4 + \lambda_1 = 0\\ \frac{\partial \mathcal{L}}{ \partial y} & = -x + 1 + \lambda_2 = 0 \\ x&=0\\ y&=1 \end{aligned} \end{equation}
hence here $x=y=0$ and $\lambda_1 = 3$ and $\lambda_2=-1$
\subsubsection*{Third Lagrangian:}
\begin{equation} \mathcal{L} = x^2 - xy + y - 4x - \lambda_1(-x) - \lambda_2 (-y) \end{equation}
with FOCs
\begin{equation} \begin{aligned} \frac{\partial \mathcal{L} }{\partial x} & = 2x - y -4 + \lambda_1 = 0\\ \frac{\partial \mathcal{L}}{ \partial y} & = -x + 1 + \lambda_2 = 0 \\ x&=0\\ y&=1 \end{aligned} \end{equation}
hence here $x=y=0$ and $\lambda_1 = 3$ and $\lambda_2=-1$
\subsubsection{Fourth Lagrangian}
\begin{equation} \mathcal{L} = x^2 - xy + y - 4x - \lambda_1(-x) - \lambda_2 (-y) \end{equation}
with FOCs
\begin{equation} \begin{aligned} \frac{\partial \mathcal{L} }{\partial x} & = 2x - y -4 + \lambda_1 = 0\\ \frac{\partial \mathcal{L}}{ \partial y} & = -x + 1 + \lambda_2 = 0 \\ x&=2\\ y&=1 \end{aligned} \end{equation}
Since $f$ is a quadratic function having a saddle point it is not convex.
In order to find the maximum of $f$ on the boundary $\partial R$ of the given rectangular domain $R$ just pull back $f$ to the four edges of $R$, i.e., consider the four auxiliary functions $$\eqalign{&\phi_0(y):=f(0,y)=y,\quad\phi_2(y)=f(2,y)=-y-4\qquad\qquad(0\leq y\leq 1),\cr &\psi_0(x):=f(x,0)=(x-2)^2-4,\quad \psi_1(x):=f(x,1)=(x-2.5)^2-5.25\qquad(0\leq x\leq2)\ .\cr}$$ Since all four of these functions are monotone on the relevant intervals we can conclude that $$\max_{(x,y)\in R}f(x,y)=\max\bigl\{f(0,0),f(2,0),f(0,1),f(2,1)\bigr\}=1\ .$$