We have studied multiple excitons (up to four electron-hole pairs) bound either by themselves or to a neutral impurity center in indirect semiconductors via a variational quantum Monte Carlo method. The ground state energies and particle density distributions in these complex systems are calculated. The calculated binding energy of each complex against the dissociation into a smaller complex plus an exciton is in fair agreement with experiment. The discrepancy is attributed to the crystal field splitting and central-cell effects present in the realistic systems.