The magnetization M and the magnetic field intensity H are conventionally used as implementational variables in the Preisach model. However, the distribution function obtained from M-H curves is strongly localized, causing difficulty during iterative computations. In this paper, a method that uses the magnetization M and the magnetic nux density B instead of M and H in the Preisach model is proposed. The Preisach distribution function and the Everett function obtained from M-B curves are compared with those of the conventional method. The former functions are much smoother than the latter, and this improves convergence and efficiency in computations because the slope of the M-B curves is lower than that of the M-H curves. A two-permanent-magnet system, as an example, is chosen and analyzed. The results are in good agreement with experiment.