We are concerned with the three-dimensional magnetostatic problem where the nonhomogeneities and the source are confined to a bounded domain. We derive a mixed formulation of this problem whose unknowns are the magnetic field, a current vector potential which is introduced as an auxiliary unknown, and a boundary unknown which results from the boundary integral method. This formulation is an improvement with respect to previous formulations proposed in the literature in the sense that it leads to an easier implementation using Nédélec's edge elements and boundary elements, and to good numerical accuracy. Some numerical results are described and compared with those obtained by using the classical formulation in a scalar potential.