The charge simulation method is an integral equation technique for solving Laplace's equation and is commonly used for potential calculations in electrostatics such as the simulation of high voltage insulation systems. The method is based on representing the potential in the exterior of a 
conductor or in dielectric regions as the superposition of a finite number of discrete fictitious  charges judiciously placed outside the region in which the field is to be computed. Although the method has enjoyed considerable popularity over the last 25 years, the computational problem it solves is inherently unstable, and much caution must be exercised in its implementation. The ill-posed nature of this problem also manifests itself when iterative techniques such as Krylov subspace methods -- which have become feasible of late due to the development of wavelet-based matrix compression techniques -- are used to solve the linear system of equations which arises in the determination of the magnitudes of the fictitious charges. The objective of this talk is to point out the ill-posed nature of the charge simulation method and to demonstrate how Krylov subspace methods can be used not only as linear solvers, but at the same time as regularization techniques to stabilize the method by terminating the iteration prematurely, before the effects of noise begin to degrade the solution.