You can probably try this quadratic function 
, where 
is a constant, and 
, 
.
, where is a constant, and , .Update 1: Another candidate function for the surface is hyperbolic function, 
.
.Update 2: If you insist on using the periodic trigonometric function on a non-oscillatory surface, then modify the fitting model to become 
so that it fits the local data.
so that it fits the local data.Update 3: Added some surface plots.
[X, Y] = meshgrid(-1:0.05:1);
%% Polynomial function
c = 1;
b = 1/2;
a = b;
Z = c - b*Y.^2 - a*X.^2; % can try higher-order polynomial functions
surf(X, Y, Z), grid on
title('Quadratic function'), xlabel x, ylabel y, zlabel f(x,y)
%% Hyperbolic function
c = 1/(1 - sech(1));
b = [1, (c*sech(1))/2];
a = b;
Z = c - b(2)*cosh(b(1)*Y) - a(2)*cosh(a(1)*X);
surf(X, Y, Z), grid on
title('Hyperbolic function'), xlabel x, ylabel y, zlabel f(x,y)
%% Trigonometric function (local)
c = 0;
b = [pi/2, 1/2];
a = b;
Z = c + b(2)*cos(b(1)*Y).^2 + a(2)*cos(a(1)*X).^2;
surf(X, Y, Z), grid on
title('Trigonometric function'), xlabel x, ylabel y, zlabel f(x,y)
%% Gaussian function (local)
c = 1/(1 - exp(1));
b = exp(1)/(2*(exp(1) - 1));
a = b;
Z = c + b*exp(-Y.^2) + a*exp(-X.^2); % can sech(x) if you like single-hump functions
surf(X, Y, Z), grid on
title('Gaussian function'), xlabel x, ylabel y, zlabel f(x,y)
