Finite difference schemes currently applied to the modeling of two-phase flows in flow networks exhibit difficulties in properly simulating certain spatial and temporal discontinuities. These discontinuities include points along the one-dimensional flow axis where density and other thermophysical properties become discontinuous or experience rapid state domain changes. A methodology for treating spatial and temporal discontinuities is presented. This methodology consists of three main features: (a) subnode time-averaged do-noring of thermodynamic properties, (b) a variable pressure-at-discontinuity staggered mesh discretization, and (c) a variable point state equation linearization. The proposed scheme is similar in form to standard semi-implicit, staggered mesh discretizations, requires little extra overhead, and results in substantially improved accuracy and code execution times. Comparisons are made with standard time and spatial discretizations, as well as with two simpler alternate methods for recognizing and tracking discontinuities. The first of these attempts is to adjust the time-step size such that the fluid discontinuity arrives at a node boundary, or a change in fluid state occurs precisely at the end of a time advancement. The second attempts to redistribute mass and energy to correct for improperly donored values when a discontinuity crosses a node boundary during a time step. Neither of these alternatives proved adequate.