One thing you can do is compute the Hessian:
[tex]\mathbf H(x,y)=\begin{bmatrix}f_{xx}&f_{xy}\\f_{yx}&f_{yy}\end{bmatrix}=\begin{bmatrix}-2&0\\0&-2\end{bmatrix}[/tex]
We have [tex]\det\mathbf H(x,y)=4>0[/tex], and [tex]f_{xx}(x,y)=-2[/tex] for all points [tex](x,y)[/tex], which means (by the second partial derivative test) that any critical point of [tex]f(x,y)[/tex] is the site of a local maximum.
The critical points occur where [tex]f_x=f_y=0[/tex]. We have
[tex]f_x=4-2x=0\implies x=2[/tex]
[tex]f_y=6-2y=0\implies y=3[/tex]
so that (2, 3) is the only critical point of [tex]f[/tex], and it happens to fall within the region [tex]D[/tex]. At this point, we have [tex]f(2,3)=20[/tex].
Meanwhile, on the boundary:
If [tex]x=0[/tex], then [tex]f(0,y)=6y-y^2=9-(y-3)^2[/tex], which has a maximum of 9 when [tex]y=3[/tex], and a minimum of 0 when [tex]y=0[/tex]. Then [tex]f(0,3)=16[/tex] and [tex]f(0,0)=7[/tex].
If [tex]x=4[/tex], then [tex]f(4,y)=7+6y-y^2=16-(y-3)^2[/tex], which has a maximum of 16 when [tex]y=3[/tex], and a minimum of 5 when [tex]y=0[/tex]. Then [tex]f(4,3)=16[/tex] (and we already know [tex]f(0,0)=7[/tex]).
If [tex]y=0[/tex], then [tex]f(x,0)=7+4x-x^2=11-(x-2)^2[/tex], which has a maximum of 11 when [tex]x=2[/tex], and a minimum of 0 when [tex]x=0[/tex] or [tex]x=4[/tex]. Then [tex]f(2,0)=11[/tex].
If [tex]y=5[/tex], then [tex]f(x,5)=12+4x-x^2=16-(x-2)^2[/tex], which has a maximum of 16 when [tex]x=2[/tex], and a minimum of 12 when [tex]x=0[/tex] or [tex]x=4[/tex]. Then [tex]f(2,5)=16[/tex].
So, over the region [tex]D[/tex], [tex]f(x,y)=4x+6y-x^2-y^2+7[/tex] has an absolute maximum of 20 at (2, 3), and an absolute minimum of 7 at (0, 0).
