A new numerical method is introduced for the solution of Poisson's equation for the electrostatic potential between arbitrarily shaped boundary surfaces that may appear in metal-molecule-metal junctions. This method is based on a straightforward procedure in which the arbitrarily shaped system is embedded in a cubic box. The embedding procedure is formulated in terms of boundary operators that can be readily implemented even for complex irregular geometries of the boundary surfaces. The solution to Poisson's equation on a cubic mesh (i.e., the inverse Laplacian operation) is used as a preconditioner, and the solution of the noncubic, more complex electrostatic problem is obtained by an error-minimization scheme that is based on a Krylov subspace expansion method. The accuracy and fast convergence of this numerical procedure are demonstrated for generic examples.