A new iterative method for computing three-dimensional steady-state magnetic fields with eddy currents is presented. By using the proposed method, the numerical computation of eddy current fields can be divided into two successive stages on flux density and eddy current calculations. The convergent field solution is then obtained iteratively. The coefficient matrices arising from the proposed method contain relatively few variables and are real. As these matrices need to be eliminated only once in the iteration procedure, the requirement upon the computer resource can be reduced substantially. The convergence of the presented iterative method is also discussed in detail. The instructions for choosing the penalty factor and relaxation factor in order to obtain the globally convergent potentials with sufficiently accurate field solutions are also given. Some sample calculations show that the new iterative method is highly computationally efficient for studying large-scale unbounded eddy-current problems in engineering.