% Grid
y=70:160;
x=0:60;
z=10:80;
[X,Y,Z]=meshgrid(x,y,z);
% Evaluate function at grid points
F=-1.69*X +.78*Y + .739*Z - .0165*X.^2 - .0106*Y.^2 + .0188*X.*Y + ...
.000097*X.^3 + .000039*Y.^3 - .000151*X.^2.*Y - .000105*X.*Y.^2 + .000203*X.*Y.*Z - 6;
% Iso-values of F
n=30; % number of evenly spaced isocontours to visualize
F_min=min(F(:));
F_max=max(F(:));
df=(F_max-F_min)/(n+2);
v=round((1:n)*df);
% Visulize xy-slices at z= 10, 30, 50, and 70
figure('color','w')
zi=[10 30 50 70];
id=[1 21 41 61]; % xy slice indices
for i=1:numel(id)
ha=subplot(1,4,i);
axis equal on
hold on
imagesc([x(1) x(end)],[y(1) y(end)],F(:,:,id(i)));
set(ha,'Xlim',[x(1) x(end)],'Ylim',[y(1) y(end)],'box','on','CLim',[F_min F_max])
[~,h]=contour(X(:,:,id(i)),Y(:,:,id(i)),F(:,:,id(i)),v,'LevelListMode','manual','ShowText','on');
set(h,'LineColor','k','LineWidth',1)
uistack(h,'top')
set(get(ha,'Title'),'String',sprintf('z = %.1f',zi(i)),'FontSize',15)
end
