Neil - the answer is definitely yes, and no. Or maybe it is the other way around.
No, it would be difficult to make gridfit work as you would like to see via internal changes, since that would be a bound constraint on a fairly massively large scale linear algebra problem. So while I could in theory use a tool like lsqnonneg or lsqlin instead of backslash for the solve, it would take forever to terminate on any problem of reasonable size.
Regardless, there is a solution that will work for you, and is quite simple. Typically when there is a bound constraint at zero because of physical issues like this it is because the problem really should be transformed to make it more linear. The logical transformation is a log transform. Thus instead of fitting Z as a function of X and Y, fit the surface in the form of log(Z) as a function of X and Y. Clearly there will no longer be any negativity issues, because to recover the surface you really want, you will simply exponentiate.
If you actually had any data points that were an exact zero, it would be best to drop them completely. Alternatively you could replace them with some small number, but that would cause the problem to be less smooth and might introduce bugs in the surface. Of course, the non-smooth region will be in areas where the log will be very negative, so after exponentiation it will be essentially zero anyway. So either way will work on those exact zeros.
I hope the log transform idea solves your problem. It is an approach I have used often with good success.

Thank you for making this code available. Is there a way to tweak it to add the additional constraint that the fitted surface should also always be positive? The X,Y,Z data I'm fitting represent spectrum (Z) vs. energy (Y) and angle (X). Therefore, negative Z is unphysical.
Thanks,
Neil

Not only a great tool, the source code and the discussions here are a great education. Thanks, John, for extending the docs and examples
It usually does everything I need. When really off-the-wall needs have come up, it's been a great base to start from. Clear and well structured code.

Hi Yair - My code snippet replaced my earlier comment. The 'aliasing' problem is that findjobj cannot distinguish two different tables placed in the same position in two different tabs.
Thanks!

Here is the test code to demonstrate the aliasing problem:
% findjobj test script
hWnd=figure;
tgrp = uitabgroup('Parent',hWnd);
tab1 = uitab('Parent',tgrp,'Title','tab1');
tab2 = uitab('Parent',tgrp,'Title','tab2');
t1=uitable('Parent',tab1);
a=findjobj(t1),
%%
t2=uitable('Parent',tab2);
%%
a=findjobj(t1),
b=findjobj(t2),
a==b,

Comment only