The 'poly2' option for fit will use a simple linear least squares solver. (I did verify this fact, as could you have done. Be VERY CAREFUL, if you edit the code to view it, as it can be a dangerous thing. Far too many people have editted code from MATLAB, and then mistakenly introduced bugs into the code, and then saved the file by mistake. Use type instead to view provided code.) As well, you could use the debugger to step down into the code to see exactly what fit does and where it goes for the fit.
Specifically, it creates the appropriate matrix problem for a quadratic model, then uses MATLAB's economy sized QR decomposition to decompose the matrix in an efficient and numerically well posed form. Why a QR? Because this also allows the computation of variances on the parameters, as will be computed later on. Essentially, R then provides a Cholesky factor for the covariance matrix. QR is a numerically stable solver for this class of linear agebra problems.
No column pivoting is performed in the QR. This is the only part that mildly surprised me, in that column pivoting could be performed to make the solve somewhat more stable. However, if your model is that close to the edge of instability that pivoting would help here, then you are trying to fit too high order a polynomial model in the first place. As such, I can accept the choice of no column pivoting done in the solve. (By way of comparison, when I wrote my own polyfitn utility, I did use a column pivoted QR.)
As I said, this is just simple linear least squares though. I'm not sure what you are asking, if it is beyond what I've said here. You can ask if you have some further question on what is done, though I won't be willing to teach a complete course on numerical linear algebra.
