NLINFIT is good.
NonLinearModel.fit is also good.
But since your problem involves fitting a surface with only two independent variables, it can be done very simply using the Curve Fitting Toolbox functions.
You can do it interactively using CFTOOL and then generate the MATLAB code automatically (recommended), or if you want to write the code out by hand yourself, you can do something along these lines:
x = rand(100,1);
z = rand(100,1);
atrue = 2.5;
btrue = 1.7;
ctrue = 1.2;
Y = atrue*(x.^btrue).*(z.^ctrue) + 0.05*randn(size(x)) ;
scatter3(x,z,Y);
hold all;
F = fittype('a*x^b*z^c','Independent',{'x' 'z'});
M = fit([x z],Y,F) % Or specify an initial guess: M = fit([x z],Y,F,'Start',[0 0 0])
plot(M)
