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.

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.