In liquid-liquid contacting equipment such as completely mixed and differential contactors, droplet population balance based modeling is now being used to describe the complex hydrodynamic behavior of the dispersed phase. For the hydrodynamics of these interacting dispersions this model accounts for droplet breakage, droplet coalescence, axial dispersion, exit and entry events. The resulting population balance equations are integro-partial differential equations (IPDE) that rarely have an analytical solution, especially when they show spatial dependency, and hence numerical solutions are sought in general. To do this, these IPDEs are projected onto a system of convective dominant partial differential equations by discretizing the droplet diameter (internal coordinate). This is accomplished by generalizing the fixed-pivot (GFP) technique of Kumar and Ramkrishna (Chem. Eng. Sci. 51 (1996a) 1311) handling any two integral properties of the population number density for continuous flow systems by treating the inlet feed distribution as a source term. Moreover, the GFD technique has the advantage of being free of repeated or double integral evaluation resulting from the weighted residual approaches such as the Galerkin's method. This allows the time-dependent breakage and coalescence functions to be easily handled without appreciable increase in the computational time. The resulting system of PDEs is spatially discretized in conservative form using a simplified first order upwind scheme as well as first- and second-order non-oscillatory central differencing schemes. This spatial discretization avoids the characteristic decomposition of the convective flux based on the approximate Riemann solvers and the operator splitting technique required by classical upwind schemes. The time variable is discretized using an implicit strongly stable approach that is formulated by careful lagging of the non-linear parts of the convective and source terms. The algorithm is tested against analytical solutions of the simplified population balance equation for a differential liquid-liquid extraction column through four case studies. In all these case studies the discrete models converge successfully to the available analytical solutions and to solutions on relatively fine grids when the analytical solution is not available. Realization of the algorithm is accomplished by comparing its predictions to experimental steady-state hydrodynamic data of a laboratory scale rotating disc contactor of 0.15 m diameter. Practically, the combined algorithm is found fast enough for the computation of the transient and steady-state hydrodynamic behavior of the continuously and spatially distributed interacting liquid liquid dispersions. (C) 2004 Elsevier Ltd. All rights reserved.