Code from the Appendix

Supporting Interval Functions

mid[X_] := (Min[X] + Max[X])/2;

Subdivision of Intervals and Rectangles

subdivide1D[X_] := 
        Interval /@ {{Min[X], mid[X]}, {mid[X], 
              Max[X]}};
subdivide2D[{X_, Y_}] := 
        Distribute[{subdivide1D[X], subdivide1D[Y]}, List];

Two-Dimensional Minimization

LowestCriticalPoint[f_,  {x_,  a_,  b_},  {y_,  c_,  d_}, 
          upperbound_,  
          tol_] := (rects = 
            N[{ {Interval[{a, b}], 
                  Interval[{c, d}]}}]; 
fcn[{xx_, 
                yy_}]  :=  
            f /. { x ->  xx, 
                y ->  yy}; 
gradf[{xx_, yy_}] := 
            Evaluate[{D[f, x], D[f, y]} /. {x -> xx,  
                  y ->  yy}]; 
{low, 
              upp} = {-∞, 
              upperbound}; 
While[(upp - low) > 
              tol, 
rects = 
              Join  @@  subdivide2D  /@  
                  rects; 
fvals = 
              fcn  /@  rects; 
upp = 
              Min[upp, Min[Max /@ fvals]]; 
pos = 
              Flatten[Position[
                  Min /@ fvals, _?(# ≤ 
                          upp  &)]]; 
rects = 
              rects[[pos]];  
            fvals = fvals[[
                pos]]; 
pos =  
              Flatten[Position[
                  Apply[And, IntervalMemberQ[gradf /@  rects, 0],  {1}],  
                  True]]; 
rects = 
              rects[[pos]]; 
            low = Min[
                fvals[[
                  pos]]]]; 
{{low, upp}, 
            Map[mid,  rects,  {2}]});

Definition of f

f[x_, y_] := (E^Sin)[50 x] + 
          Sin[60 (E^y)] + Sin[70 Sin[x]] + Sin[Sin[80 y]] - 
          Sin[10 (x + y)] +  ((x^2) + (y^2))/4;
f[{x_, y_}] := f[x, y];

Session p. 85

LowestCriticalPoint[f[x, y], {x, -1, 1}, {y, -1, 1}, -3.24, 
      (10^(-9))]