Hey Grega,
I'm glad you asked this question, because I now know 2 different solutions. The first solution is the manual solution that I worked out myself, which I won't give here. The second is from a Mathworks workshop I went to, and looks like this:
X = x2fx(predictors,'quadratic');
y = responses;
[Q,R] = qr(X,0);
beta = R\(Q'*y);
yhat = X*beta;
residuals = y - yhat;
where predictors is a p by m vector or row matrix of data and responses is a p by 1 vector. This is all given in the code for regstats. I should have looked there first it seems.
Best,
Sam