Kriging is used extensively as a metamodel in multidisciplinary design optimization. The correlation matrix used in Kriging metamodeling frequently becomes ill-conditioned. Therefore different numerical methods used to solve the Kriging equations affect the search for the optimum Kriging parameters and the ability of the Kriging surface to accurately interpolate known data points. We illustrate this by firstly computing the inverse of the correlation matrix in the Kriging equations, and secondly by solving the systems of equations using decomposition and back substitution, thereby avoiding the inversion of the correlation matrix. Our results clearly show that by decomposing and back substituting, the interpolation accuracy is maintained for significantly higher condition numbers. We then show that computing the natural logarithm of the determinant using additive calculations as opposed to multiplicative calculations significantly reduces numerical underflow errors encountered when searching for the optimum Kriging parameters. Although the effect of decomposition and back substitution are known, and the underflow difficulties when computing the natural logarithm of the determinant of the correlation matrix has been mentioned in passing in Kriging literature, this work clearly quantifies and reinforces these methods, hopefully for the benefit of researchers entering the field.