%  You need to download   www.glue.umd.edu/~tvp/241/ezcontourc.m

syms x y
f = -x^4 - y^4 - x^2*y^2 + 6*x^2 + 6*y^2

fx = diff(f,x)                           % find first partial derivatives of f
fy = diff(f,y)

fxx = diff(fx,x); fxy = diff(fx,y);      % find second partial derivatives of f
fyy = diff(fy,y)
D = fxx*fyy - fxy^2                      % find D

[xc,yc] = solve(fx,fy,x,y); [xc,yc]      % solve fx=0, fy=0 for x,y
nc = length(xc);                         % Matlab finds only 7 of 9 solutions

for i=1:nc                               % for each critical point
  fprintf('Critical point %i = \n',i)
  disp([xc(i),yc(i)])
  fxx_value = subs(fxx,{x,y},{xc(i),yc(i)})  % evaluate fxx at critical pt.
  D_value =   subs(D,{x,y},{xc(i),yc(i)})    % evaluate D   at critical pt.
end

figure(1); ezsurf(f,[-2,2,-2,2])         % surface plot 
figure(2); ezcontourc(f,[-2,2,-2,2],15)  % contour plot with 15 contours
axis equal; axis tight; colorbar
hold on

plot(double(xc),double(yc),'kx')         % mark critical points with black "x"

hold off



