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.