  The equations of state for symmetric nuclear matter and pure neutron matter
are investigated with the tensor-optimized Fermi Sphere method (TOFS) up to the
density $\rho=0.5$~fm$^{-3}$. This method is based on a linked-cluster
expansion theorem, and the energy per particle of nuclear matter ($E/A$) is
calculated variationally with respect to the correlated nuclear matter wave
function. We can study the density dependence of the many-body terms arising
from the operator products, which contribute to $E/A$. In order to clarify the
relation between the many-body effects and short-range correlation, we take the
spin-isospin dependent central {\it NN} interaction with a few GeV repulsion in
the inner region. The EOS obtained by the TOFS method is reasonably reproduced,
compared with other \textit{ab initio} many-body methods. We found that the
many-body terms from the 2-body to 6-body ones) give sizable effects on $E/A$
at higher density, and they play an important role in nuclear matter.
