You will NEED to use a spline based method of interpolation here, specifically spline or cubic. NOT makima, or pchip. Those interpolants will not generate functions with an extremum between the data points. And of course you certainly cannot use the default linear interpolant, as that will NEVER generate a solution that is not at a data point.
[X,Y] = ndgrid(0:0.5:6);
Z = cos(X+Y).*sin(X-Y);
G = griddedInterpolant(X,Y,Z,'spline');
fsurf(@(x,y) G(x,y),[0 6 0 6])
We can clearly see several peaks. Of course, since I used a trig function here, they will lie at transcendental locations, so never exactly at a grid point.
[xymin,fval] = fminsearch(@(xy) -G(xy(1),xy(2)),[1,1])
To within the tolerances used by fminsearch, this should be a solution.
xymin/pi
And that would clearly be the point at [pi/4,3*pi/4]. If you want a better solution you need a finer grid spacing, as the spline has limited approximative capabilities with that course grid.
