This snippet, is in reality the initial prototype I used for building my polynomial least squares class/module in C and Fortran, respectively. As you'll see, it matches perfectly the outputs for Matlab's polyval() and lu() functions. Although, is a very minimalist implementation of the LU decomposition method, it sometimes beats the polyval() function in speed. (Not sure why?!)

To understand how to use it, check out the example : FitPolynomialToNoisyData.m
( In it, I compare this implementation to Matlab's traditional tools/approaches )

Check it out! ;D