f := (x,y) -> 2*exp(-(x^2+y^2)/2)/3 + exp(-((x-2)^2+(y+2)^2)/2)/3: plot3d(f(x,y),x=-4..4,y=-4..4,axes=frame); Globalmax := f(0,0); localmax := f(2,-2); evalf(Globalmax), evalf(localmax); fsolve({diff(f(x,y),x),diff(f(x,y),y)},{x,y}); fsolve({diff(f(x,y),x),diff(f(x,y),y)},{x,y},{x=-0.5..0.5,y=-0.5..0.5}); g := (x,y) -> x^2 - y^2: gg := grad(g(x,y),[x,y]); plot3d(g(x,y),x=-1..1,y=-1..1,axes=frame);