To be learned in school, the maxima and minima of a function are located, where the first derivative is zero.
So one has in two dimensions to task to calculate the gradient and the Hessian Matrix.
Your are lucky that this is a Mathematica help example:
FindMaximum[Sin[x] Sin[2 y], {x, y},Gradient -> {Cos[x] Sin[2 y], 2 Cos[2 y] Sin[x]}, Method -> {"Newton", Hessian -> {{-Sin[x] Sin[2 y],
2 Cos[x] Cos[2 y]}, {2 Cos[x] Cos[2 y], -4 Sin[x] Sin[2 y]}}}]
On the page for FindMaximum
.
Smallest positive solution is {x -> π/2, y -> π/4}
The example from Mathematica shows, one has to think further on ones own. Both factors of the given function are trigonometric and therefore periodical. One factor has the frequency 1, the other two. So all one has to do after the work You have already done is to add this periodical to the solution. And then restrict them to the region given.
So there is a grid of solutions starting in the positive quadrant with the solution You have calculated. Add 2 Pi or Pi to the component values of the solution and compare to 0 or 4 will to the job.
There is a speciality in the discussion not covered with periodicity.
There is a relative maximum. That You already saw on the border of the given rectangle. So to complete the job restrict the two dimensional function on the border lines.
x==0 is constant, so no maximum there.
y==0 is constant, so no maximum there.
x==4 is Sin[4]Sin[2y]=-0.756802 Sin[2y]
starts with 0 and gets negative for positive y. The maximum is smaller than y=4, so {{4,0},{4,3π/2}}
are maxima in the given region.
y==4 is Sin[x]Sin[8]=0.989358 Sin[x]
starts positive and rises for positive x.
That is the complete job and does contain the mathematical full and accepted procedure for this task.
(i) Find the minima/maxima in the open given region.
(ii) Find minima/maxima on the border of the region.
Sorry that Solve is not much better. It is just equal and saves You to think like a Mathematician.
My methodology confirms the steps already done by Bob Halon, but presents the real world maths to be done and requested by staff for example in exams. No searching fo r a good starting point is needed, just straight forward math.
The points with maxima in the region {{0,4},{0,4}}
are
{{x -> π/2, y -> π/4}, {x -> π/2, y -> (5 π)/4},{4,0},{4,3π/2}}.

The methods choosen behind the function FindMaximum
are state of the art or really fast methods. Wolfram Inc. does not offer too much information about what the function really are. If You once are in the methodology of calculating maxima and minima in Mathematical Numerics there will be only a few worth to be choosen. A nice start is From Curve Fitting to Machine Learning: An Illustrative Guide to Scientific ...
from Achim Zielesny. For further insight You have to study the references therein.
FindMaximum
with initial points to search, because it finds local maxima. $\endgroup$