See if this does what you want:
dsnmtx = [ones(size(y',1),1) x1' x2']; % Design Matrix
B = dsnmtx\y'; % Estimate Parameters
yfit = [ones(size(y',1),1) x1' x2']*B; % Estimated ‘y’
figure(1)
stem3(x1', x2', y')
hold on
plot3(x1', x2', yfit)
hold off
grid on