Continuously differentiable radial basis functions (C ∞ -RBFs), while being theoretically exponentially convergent are considered impractical computationally because the coefficient matrices are full and can become very illconditioned. Similarly, the Hilbert and Vandermonde have full matrices and become ill-conditioned. The difference between a coefficient matrix generated by C ∞ -RBFs for partial differential or integral equations and Hilbert and Vandermonde systems is that C ∞ -RBFs are very sensitive to small changes in the adjustable parameters. These parameters affect the condition number and solution accuracy. The error terrain has many local and global maxima and minima. To find stable and accurate numerical solutions for full linear equation systems, this study proposes a hybrid combination of block Gaussian elimination (BGE) combined with arbitrary precision arithmetic (APA) to minimize the accumulation of rounding errors. In the future, this algorithm can execute faster using preconditioners and implemented on massively parallel computers.