In multiphase 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 derivation is relatively simple and mainly
involves centered finite difference formulas, with less reliance on the Taylor
series expansions of typical immersed interface method derivations. 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.