I don't have the Stats Toolbox, but according to the online documentation, you can do exactly what you've shown, assuming x,y,z are all column vectors.
As for the initial parameter values, that's always a matter of ingenuity, and insight into the particular model function that you have. Often, you would look for some semi-accurate approximation to the model function, that is easier to fit, and use that simplified fit as your initial guess.
I cannot use "lsqcurvefit" as it does not give the fitting error.
That shouldn't be a major barrier, in any case. LSQCURVEFIT does return the residual and Jacobian from which your other goodness of fit measures can be easily calculated.
