The normal equations for least squares are
=MMULT(MINVERSE(MMULT(TRANSPOSE(xMat),xMat))),MMULT(TRANSPOSE(xMat),yVec))
Which will work for up to 52 predictors if your X matrix is orthogonal.
This "method of calculation ... breaks down" for non-orthogonal X
matrices in the sense that the the useful information often ends up well
beyond the accuracy of IEEE double precision. A numerically better
approach is to do a singular value decomposition of the X matrix, so
that much of the numerical junk in (X'X)^-1 X'y can be canceled
analytically instead of needing far more precision than is available.
My vague recollection is that MATLAB also forms the normal equations,
and hence will have the same numerical instability as the excel formula
above. R (which Harlan suggested) is free
http://www.r-project.org/
and uses a more stable algorithm akin to the approach that I outlined above.
Jerry