A high-order difference scheme is derived from the dilatational wave equations of two-phase media, and the corresponding absorbing boundary conditions and stability conditions are also provided, based on which a high-order finite difference numerical modeling algorithm for dilatational wave equations of two-phase isotropy media is realized. It is shown that this method can greatly improve the precision on condition that the calculation time is a bit extended. This method can be applied to both pre-stack and post-stack numerical modeling.