In the formulation with two magnetic vector potentials, the inaccuracy of the Biot-Savart integration for the source term causes the divergence of the conjugate gradient method. The ill-posed equations are regularized by subtracting a rotational field from the numerically integrated field. The convergence is drastically improved. Even when the current continuity is not satisfied strictly, we can get reasonable results. Also, the response to a small disturbance can be calculated by the regularization, even, when the modeling of source currents. is fairly rough.