  We present a novel approach for computing reduced density matrices for
superpositions of eigenstates of a Bethe-ansatz solvable model by direct
integration of the wave function in coordinate representation. A diagrammatic
approach is developed to keep track of relevant terms and identify symmetries,
which helps to reduce the number of terms that have to be evaluated
numerically. As a first application we compute with modest numerical resources
the single-particle density matrix and its eigenvalues including the condensate
fraction for a quantum bright soliton with up to $N=10$ bosons. The latter are
constructed as superpositions of string-type Bethe-ansatz eigenstates of
nonrelativistic bosons in one spatial dimension with attractive contact
interaction. Upon delocalising the superposition in momentum space we find that
the condensate fraction reaches maximum values larger than 97\% in the range of
particles studied. The presented approach is suitable for studying
time-dependent problems and generalises to higher-order correlation functions.
