There are several ways one could do this, and YES, it is possible. There are some issues one must resolve of course.
The simplest solution is to use a tool like my gridfit. That tool is simply a piecewise linear regression spline, defined in terms of the function values at each node in the grid. So were you to build a 100x100 grid, there are 10000 unknowns to estimate using a regression scheme. The constraint that dz/dx>=0 everywhere is simple to apply here of course, since one simply forms one inequality constraint equation for each pair of neighboring x-connected nodes. So there would be 99*100 such linear inequality constraints. The resulting system might be rather large for LSQLIN to solve I fear. However, you could probably solve a 20x20 problem easily enough, so only 400 unknowns, and 380 linear inequality constraints.
As I said, the above is the simple solution, although it might take some CPU time depending on how large of a grid you chose to populate.
Next, we can consider a 3rd order model in x and y. I'll assume:
z(x,y) = a00 + a10*x + a01*y + a20*x^2 + a11*x*y + a02*y^2 + a30*x^3 + a21*x^2*y + a12*x*y^2 + a03*y^3
There are 10 coefficients here. The problem, if we wish to constrain the slope, dz/dx to be positive everywhere inside the bounds, is the constraint becomes a nonlinear one. Worse, since it must be true everywhere in that domain, it gets a bit nasty. A simple solution is to define a set of NECESSARY constraints, for positivity. That is, at each of a large set of points in the (x,y) domain of interest, we will set a linear constraint such that
dz/dx >= 0
Once a point (x_k,y_k) is chosen, that derivative becomes a linear inequality constraint in the parameters.
If we choose a large enough set of points for that set of linear inequality constraints, then you can see that there MAY be some small region in the domain of interest where the constraint may SLIGHTLY fail, but only by a small amount if our test set is chosen to be large enough. The result is a linear least squares problem, with 10 unknowns and as many linear inequality constraints as you choose to generate. The more constraints you choose, then better the result.
Again, the above scenario is a NECESSARY, but NOT SUFFICIENT constraint on slope positivity. It should get you pretty close however.
Another approach can be found from the work of Fritsch and Carlson. (This is the idea that originally spawned PCHIP by the way.) Here we can form a SUFFICIENT constraint on the linear system for the coefficients of the cubic model. The problem is that the constraint system would be nonlinear but still convex, so the trick is to convert that nonlinear constraint into a set of LINEAR inequalities that ensure the solution has the desired slope.
The difference between the last two options (necessary versus sufficient) is a subtle one. he former may result in a solution that may have some small areas of slightly negative slope. The latter may possibly result in a slightly too conservative fit, but it always ensures slope positivity.
Of course, there are no codes written that I know of to directly solve any of these variations, but it is not terribly difficult. In fact, option2 is almost trivial to cobble together, and option 1 not too much harder. The downside of the polynomial model options is they are simply polynomial models, and not terribly high order ones at that. And since low order polynomials tend to be somewhat poor in their abilities to fit, there will be limits here too.
My personal suggestion is to throw together option 2. See how well it does, if there are regions of negative slope encountered. (Note that this will not generate a GLOBALLY positive slope, only a slope that is everywhere positive in the region of interest.) If the third order model was deemed insufficient after seeing the fit quality and lack of fit, then I'd suggest option 1 for a reasonably small grid, perhaps 20x20, maybe a bit larger. See how fine you could go before causing LSQLIN to have a hissy fit.
