An advanced numerical scheme based on the finite difference method is presented for the full solution of a steady, incompressible, viscous, and electrically conducting flow with relatively large Hartmann numbers and interaction parameters. The solution method solves the coupled Navier-Stokes and Maxwell's equations at a low magnetic Reynolds number, and focuses on predicting key performance parameters for the design of fusion liquid metal blankets. The computational algorithm includes the effects of advection and diffusion, and is intended as a full solution, which has the potential capability of treating unsteady flow and applicability to heat and mass transfer in MHD flows.The present method uses primitive variables (velocity and pressure) for hydrodynamic variables and the electric potential, and employs a finite volume approach with a staggered grid system for fast convergence and physical accuracy. A pressure equation of the Poisson type is used for pressure-velocity-potential coupling.Results for two-dimensional channel flows with a nonuniform magnetic field at a Hartmann number and an interaction parameter of about 103 are obtained. The extension of the present numerical scheme to a three-dimensional MHD flow has been carried out and some preliminary resutls of a straight duct flow with a rectangular cross-section are presented.