The difference between min(x1) and max(x1) is not 100, thus the vector x1fit generated via colon() is just min(x1). Consequently, YFIT is not a 2x2 matrix as expected by mesh() (or surf() for that matter)
load('data.mat')
x1 = p21;
x2 = vp;
y = ucs;
X = [ones(size(x1)) x1 x2 x1.*x2];
b = regress(y,X)
scatter3(x1,x2,y,'filled')
hold on
%Values for verification
min(x1)
max(x1)
%Generate 50 points in between min and max values using linspace()
x1fit = linspace(min(x1), max(x1), 50);
x2fit = linspace(min(x2), max(x2), 50);
[X1FIT,X2FIT] = meshgrid(x1fit,x2fit);
YFIT = b(1) + b(2)*X1FIT + b(3)*X2FIT + b(4)*X1FIT.*X2FIT;
mesh(X1FIT,X2FIT,YFIT)
xlabel('P21 [m/m2]')
ylabel('Vp [m/s]')
zlabel('UCS [MPa]')
view(50,10)
hold off
[~,~,~,~,stats] = regress(y,X)
