<p>Accurate electrostatic field computation near geometrically complex boundaries (corners, high-curvature regions) is critical for microelectronics, high-voltage engineering, and materials science. Traditional numerical methods, such as the Finite Element Method (FEM) and Boundary Element Method (BEM), face challenges in accuracy and efficiency for such problems. This paper introduces a novel, robust multi-stage conformal mapping framework for high-fidelity electrostatic analysis. The framework systematically transforms complex physical domains and quasi-3D symmetric structures to the standard upper half-plane via cascaded conformal mappings, enabling analytical solutions via the method of images. Our core contribution lies in rigorously analyzing and resolving severe numerical instabilities inherent in the boundary integral formulation. We identify rotational non-uniqueness as a primary source of ill-conditioning (condition number ~ 10¹⁶) and propose a “point-fixing constraint” reformulation that reduces it to near-ideal levels (~ 2). Furthermore, we develop a randomized Nyström preconditioning technique to effectively suppress ill-conditioning induced by geometric singularities (corners). Numerical experiments demonstrate the framework’s superiority. For challenging 2D and 3D problems featuring complex boundaries with corners and high-curvature, our method achieves up to a 33-fold computational speedup compared to standard BEM, while maintaining high accuracy (relative L2 error norm &lt; 10⁻²). This significant efficiency gain, stemming from the semi-analytical nature of the approach and drastic reduction in degrees of freedom, establishes the framework as a powerful tool for rapid design iteration and simulation-based optimization in engineering applications.</p>