  Elliptic interface problems whose solutions are $C^0$ continuous have been
well studied over the past two decades. The well-known numerical methods
include the strongly stable generalized finite element method (SGFEM) and
immersed FEM (IFEM). In this paper, we study numerically a larger class of
elliptic interface problems where their solutions are discontinuous. A direct
application of these existing methods fails immediately as the approximate
solution is in a larger space that covers discontinuous functions. We propose a
class of high-order enriched unfitted FEMs to solve these problems with
implicit or Robin-type interface jump conditions. We design new enrichment
functions that capture the imposed discontinuity of the solution while keeping
the condition number from fast growth. A linear enriched method in 1D was
recently developed using one enrichment function and we generalized it to an
arbitrary degree using two simple discontinuous one-sided enrichment functions.
The natural tensor product extension to the 2D case is demonstrated. Optimal
order convergence in the $L^2$ and broken $H^1$-norms are established. We also
establish superconvergence at all discretization nodes (including exact nodal
values in special cases). Numerical examples are provided to confirm the
theory. Finally, to prove the efficiency of the method for practical problems,
the enriched linear, quadratic, and cubic elements are applied to a multi-layer
wall model for drug-eluting stents in which zero-flux jump conditions and
implicit concentration interface conditions are both present.
