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);
f=3.*x.*exp(y)-x.^3-exp(3.*y);
figure(1)
surf(x,y,f)
xlabel x
ylabel y
view(136,3)
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); % Derivative
figure(1);
contour(x,y,fa1, [0; 0]);
hold on;
contour(x,y,fa2, [0; 0]); % maximum where df/dx = df/dy = 0
grid on % we find one such point at (1,0]
hold off
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])
