  In multi-phase fluid flow, fluid-structure interaction, and other
applications, partial differential equations (PDEs) often arise with
discontinuous coefficients and singular sources (e.g., Dirac delta functions).
These complexities arise due to changes in material properties at an immersed
interface or embedded boundary, which may have an irregular shape.
Consequently, the solution and its gradient can be discontinuous, and numerical
methods can be difficult to design. Here a new method is presented and
analyzed, using a simple formulation of one-dimensional finite differences on a
Cartesian grid, allowing for a relatively easy setup for one-, two-, or
three-dimensional problems. The method preserves a sharp interface with
discontinuous solutions, obtained from a small number of iterations
(approximately five) of solving a symmetric linear system with updates to the
right- hand side. Second-order accuracy is rigorously proven in one spatial
dimension and demonstrated through numerical examples in two and three spatial
dimensions. The method is tested here on the variable-coefficient Poisson
equation, and it could be extended for use on time-dependent problems of heat
transfer, fluid dynamics, or other applications.
