The main problem with your plot is simply that the large values of f for y > 0.5 dwarf the variations at lower values. Changing the plotting range solves this problem, I also rotate the axes to make the graph more similar to the original.
[x,y] = meshgrid(-1.5:.1:1.5,-2:0.1:0.5);
The critical points are where both partial derivatives are zero. You can visually find one solution by plotting the zero contours for bot in the same graph:
[fa1,fa2] = gradient(f, 0.2, 0.2);
contour(x,y,fa1, [0; 0]);
contour(x,y,fa2, [0; 0]);
You can also find analytical expressions for the partial derivatives and solve using fsolve from the optimization toolbox. Here I use an anonymous function J. z(1) is x and z(2) is y.
J = @(z) [3*exp(z(2))-3*z(1)^2;3*z(1)*exp(z(2))-3*exp(3*z(2))]
zc = fsolve(J,[2;0])