A numerical method was developed for predicting current density distribution in electrochemical systems of several species at steady-state. The fundamental transport equation consisted a partial differential equation (PDE) involving linear terms of diffusion and laminar convection, and nonlinear terms of ionic migration. The boundary conditions (BCs) consisted also PDEs including flux conditions and involving nonlinear terms associated with exponential kinetics of heterogeneous electrode reactions. The method of finite-difference (FD) was used to approximate solution of the global PDE system in two dimensions using a nonuniform rectangular mesh of fine size near electrode edges where boundary layers and steeper gradients were developed. Diffusion and migration terms were discretized by central finite differences, while the convective term was discretized using an upwind differencing scheme in conjunction with introduction of a 'minimal artificial viscosity' to ensure convergence and prevent anomalies in solution at high velocities. The false boundaries technique treated flux BCs. The method was tested by application to a model and an experimental system placed between two equal size parallel-plate electrodes fixed in a channel's walls under laminar parabolic flow. Numerical results for local cathodic current density were compared to experimental measurements and analytical predictions for limiting conditions. Excellent agreement was demonstrated thus verifying applicability of this FD method to real situations. (C) 2003 Elsevier Ltd. All rights reserved.