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.