Your problem is the algorithm by which you're calculating the coefficients and the variance-covariance matrix.

While it's (usually) algebraically fine, it's not numerically stable. 

You will need to use a reasonably stable algorithm to do the calculations -- they way you're doing it really can be very inaccurate; it's quite possible for *some* values to be close, and for others to be nowhere near the correct least squares solution.

A common choice in a lot of regression routines is via QR decomposition of the X matrix; other choices make more sense in particular circumstances.