The method of choice for describing attractive quantum systems is Hartree−Fock−Bogoliubov(HFB) theory. This is a nonlinear model which allows for the description ofpairing effects, the main explanation for the superconductivity ofcertain materials at very low temperature. This paper is the first study ofHartree−Fock−Bogoliubov theory from the point of view of numerical analysis. We start bydiscussing its proper discretization and then analyze the convergence of the simple fixedpoint (Roothaan) algorithm. Following works by Cancès, Le Bris and Levitt for electrons inatoms and molecules, we show that this algorithm either converges to a solution of theequation, or oscillates between two states, none of them being solution to the HFBequations. We also adapt the Optimal Damping Algorithm of Cancès and Le Bris to the HFBsetting and we analyze it. The last part of the paper is devoted to numerical experiments.We consider a purely gravitational system and numerically discover that pairing alwaysoccurs. We then examine a simplified model for nucleons, with an effective interactionsimilar to what is often used in nuclear physics. In both cases we discuss the importanceof using a damping algorithm.