Numerical simulation of two-phase flow with imbedded interface of discontinuity is challenging in two major aspects. First, the interface can undergo change, merge and breakup during the course of the simulation. Examples of current successful methods in modeling flow with interface are, among others, the volume-of-fluid method, fronttracking method, and level-set method. Second, flow variables and their derivatives can be discontinuous across the interface. This discontinuity poses severe limitation on the accuracy of commonly used numerical methods. Current available methods in treating the interface discontinuity, such as the immersed boundary method and the ghost fluid method, are mostly first order accurate. Though the immersed interface method of LeVeque and Li (1994) can be globally second order, it is often difficult to apply to complex multi-dimensional flow. This paper presents a new high-order immersed interface method for elliptic equations with imbedded interface of discontinuity. The new method can be arbitrarily high-order accurate, and it can be easily applied to practical two-phase flow because only the physical jump conditions for variables and their first derivatives are needed in the finite difference formulas. In addition, the new interface difference formulas are expressed in a general explicit form so that they can be applied to different multi-dimensional problems without any modification. The new interface algorithms of up to accuracy have been tested for one and two-dimensional elliptic equations with imbedded interface. The extension to practical two-phase flow applications will be presented in a future paper.