データ、コードを一緒に入れていただくと、回答者が取り組みやすくなります。
こんな感じで入れてみては？
x = randn(100,1);
y = randn(100,1);
z = randn(100,1);
scatter3(x,y,z);
データをフィットさせます
Fobj = fit([x,y],z,'nearestinterp'); % 最近傍法による曲面近似, 色々とやり方が有ります
% 新しいグリッド
[Xi,Yi] = meshgrid(linspace(min(x),max(y),50),linspace(min(y),max(y),50));
Zi = Fobj(Xi,Yi); % 推定値
ここまでくれば積分も出来る様に思えます:
surf(Xi,Yi,Zi);
