For the coupled system with moving boundary values of multilayer dynamics of fluids in porous media, a characteristic finite difference fractional step scheme applicable to the parallel arithmetic is put forward. Some techniques, such as the change of regions, the calculus of variations, the piecewise threefold quadratic interpolation, the multiplicative commutation rule of difference operators, the decomposition of high order difference operators, and the prior estimates, are adopted. The optimal order estimates in the l2 norm are derived to determine the error in the approximate solution. This numerical method has been successfully used to simulate the flow of migration-accumulation of the multilayer percolation coupled system. Some numerical results are well illustrated in this paper.