Custom 2D surface fitting issue

I am trying to develop some code to fit a surface. However, the values I recieve back from Matlab's fit function are horrific and non-sensical. Despite it being, in my view, an easy thing to approximate.

The surface is a summation of sine and cosine modes (here just 1 of each) and yet the surface cannot see to be recovered by the fit function.

%% 2D surface generation
clear all ; clc; close all
N=128;
L=10;
M=1; %number of modes
kmax=N/2;
kx=(randi(kmax,1,M)-kmax/2);%random wave numbers
ky=(randi(kmax,1,M)-kmax/2);%random wave numbers
x=L*(0:N-1)/N;
y=x';
V=zeros(N,N);
test = zeros(N,M);
A=rand(M,M)*2-1; %random amplitudes
B=rand(M,M)*2-1;
for i=1:M
    for j=1:M
      V=V+A(i,j)*sin((2*pi/L)*kx(i)*x+(2*pi/L)*ky(j)*y)+B(i,j)*cos((2*pi/L)*kx(i)*x+(2*pi/L)*ky(j)*y);
    end
end
figure(1)
imagesc(x,y,V)

%% Fitting
strval = [];
coeffval = {};
fac = ['(2*pi/',num2str(L),')*'];
for i=1:M
    for j=1:M
    charMi = num2str(i);
    charMj = num2str(j);
    if i==1 & j==1
    form = ['A_',charMi,charMj,'*sin(',fac,'x*','wx_',charMi,'+ ',fac,'y*','wy_',charMj,') + B_',charMi,charMj,'*cos(',fac,'x*','wx_',charMi,'+ ',fac,'y*','wy_',charMj,')'];
    else
    form = [' + A_',charMi,charMj,'*sin(',fac,'x*','wx_',charMi,'+ ',fac,'y*','wy_',charMj,') + B_',charMi,charMj,'*cos(',fac,'x*','wx_',charMi,'+ ',fac,'y*','wy_',charMj,')'];
    end
    strval = strcat(strval,form);
    coeffs = {['A_',charMi], ['A_',charMj], ['B_',charMi], ['B_',charMj], ['wx_',charMi], ['wx_',charMj], ['wy_',charMi], ['wy_',charMj]};
    coeffval = [coeffval,coeffs];
    end
end
optsL = [-ones(2*M*M,1)',zeros(2*M,1)'];
optsH = [ones(2*M*M,1)',kmax.*ones(2*M,1)'];
optsS = [zeros(2*M*M,1)',(kmax/2).*ones(2*M,1)'];

%% Fit: 'untitled fit 1'.
[xData, yData, zData] = prepareSurfaceData( x, y, V );

% Set up fittype and options.
ft = fittype( strval, 'independent', {'x', 'y'}, 'dependent', 'z' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'off';
opts.Lower = optsL;
opts.MaxIter = 5000;
opts.MaxFunEvals = 5000;
opts.Robust = 'LAR';
opts.StartPoint = optsS;
opts.Upper = optsH;
opts.Normal = 'on';

% Fit model to data.
[fitresult, gof] = fit( [xData, yData], zData, ft, opts );

% Plot fit with data.
figure(2, 'Name', 'untitled fit 1' );
h = plot( fitresult, [xData, yData], zData );
legend( h, 'untitled fit 1', 'V vs. x, y', 'Location', 'NorthEast' );
% Label axes
xlabel x
ylabel y
zlabel V
grid on
view( 9.7, 45.2 );
fitresult

Any help would be greatly appreciated.

Thanks in advance :)