This paper deals with a fully implicit finite difference scheme for the numerical solution of the Lagrangian form of the porous media fractional flow equation. It covers the theoretical background, description of mathematical formulations, basic assumptions in the models, normalization and transformation between Lagrangian and Eulerian formulations, specification of boundary conditions, discretization and solution method. Some numerical examples are described and the results compared favourably with co-current immiscible displacement data. The solution presented is more efficient than the finite element solution method in terms of requiring less matrix computation effort and it is much more stable than the explicit finite difference solution methods that have been presented in the literature. Due to the ease of formulation and use, the simulation algorithm presented can be used to formulate laboratory numerical simulators that can be used routinely for co-current flow numerical studies.