This works:
z=[0.0892 .158 0.890 1.26 2.47 3.64 6.71 11.6 18.1 27.7 42.0 65.4 97.5 141.0 196.0 283.0 440.0 661.0 863];
x=[0.05 0.628 0.0791 0.0995 0.125 0.158 0.199 0.25 0.315 0.396 0.5 0.628 0.79 0.995 1.25 1.58 1.99 2.5 3.15];
y=[1.42 1.78 2.25 2.83 3.56 4.49 5.65 7.08 8.79 11.01 13.87 17.46 21.8 27.16 33.75 42.34 53.53 67.03 84.11];
xy = [x; y];
fz= @(m,xy) 2.*m.*xy(1,:).*xy(2,:);
CF = @(m) sum((z - fz(m,xy)).^2);
m = fminsearch(CF, 1);
ze = fz(m, xy);
figure(1)
plot3(x, y, z, '+b')
hold on
plot3(x, y, ze, '-r')
hold off
grid on
view([30 40])
It estimates: m = 1.794 and seems visually to give a good fit.
