This can be accomplished using functions from the Curve Fitting Toolbox and the Global Optimization Toolbox.
The idea is to initially make a guess of the spline surface using random coefficients and the SPMAK function. Using this spline, the directional derivatives can be found with the given X and Y co-ordinates and the datum Z co-ordinate using the FNVAL and the FNDIR function. Then, using LSQCURVEFIT, we can iteratively arrive at a new guess of the spline with derivatives closer to the actual derivatives.
The sample code attached demonstrates this.
Best Answer