The ground state of a general pairing Hamiltonian for a finite nuclear system is constructed as a product of collective, real, distinct pairs. These are determined sequentially via an iterative variational procedure that resorts to diagonalizations of the Hamiltonian in restricted model spaces. Different applications of the method are provided that include comparisons with exact and projected BCS results. The quantities that are examined are correlation energies, occupation numbers and pair transfer matrix elements. In a first application within the picket-fence model, the method is seen to generate the exact ground state for pairing strengths confined in a given range. Further applications of the method concern pairing in spherically symmetric mean fields and include simple exactly solvable models as well as some realistic calculations for middle-shell Sn isotopes. In the latter applications, two different ways of defining the pairs are examined: either with J=0 or with no well-defined angular momentum. The second choice reveals to be more effective leading, under some circumstances, to solutions that are basically exact.