  We present the results of extensive computer simulations performed on
solutions of monodisperse charged rod-like polyelectrolytes in the presence of
trivalent counterions. To overcome energy barriers we used a combination of
parallel tempering and hybrid Monte Carlo techniques. Our results show that for
small values of the electrostatic interaction the solution mostly consists of
dispersed single rods. The potential of mean force between the polyelectrolyte
monomers yields an attractive interaction at short distances. For a range of
larger values of the Bjerrum length, we find finite size polyelectrolyte
bundles at thermodynamic equilibrium. Further increase of the Bjerrum length
eventually leads to phase separation and precipitation. We discuss the origin
of the observed thermodynamic stability of the finite size aggregates.
