  We have studied the properties of a classical $N_S$-body system coupled to a
bath containing $N_B$-body harmonic oscillators, employing an $(N_S+N_B)$ model
which is different from most of the existing models with $N_S=1$. We have
performed simulations for $N_S$-oscillator systems, solving $2(N_S+N_B)$
first-order differential equations with $N_S \simeq 1 - 10$ and $N_B \simeq 10
- 1000$, in order to calculate the time-dependent energy exchange between the
system and the bath. The calculated energy in the system rapidly changes while
its envelope has a much slower time dependence. Detailed calculations of the
stationary energy distribution of the system $f_S(u)$ ($u$: an energy per
particle in the system) have shown that its properties are mainly determined by
$N_S$ but weakly depend on $N_B$. The calculated $f_S(u)$ is analyzed with the
use of the $\Gamma$ and $q$-$\Gamma$ distributions: the latter is derived with
the superstatistical approach (SSA) and microcanonical approach (MCA) to the
nonextensive statistics, where $q$ stands for the entropic index. Based on
analyses of our simulation results, a critical comparison is made between the
SSA and MCA. Simulations have been performed also for the $N_S$-body ideal-gas
system. The effect of the coupling between oscillators in the bath has been
examined by additional ($N_S+N_B$) models which include baths consisting of
coupled linear chains with periodic and fixed-end boundary conditions.
