  In this paper, an approach for neutral Higgs bosons search is described based
on 2HDM type-I at electron-positron linear colliders operating at $ \sqrt{s}=1$
TeV. The beam is assumed to be unpolarized and fast detector simulation is
included. The signal process produces a fully hadronic final state through $
e^{+}e^{-}\rightarrow AH\rightarrow b\bar{b}b\bar{b} $ where both CP-even and
CP-odd Higgs bosons ($H$ and $A$) are assumed to decay to a pair of $b$-jets.
Several benchmark scenarios are introduced as the baseline for the analysis
taking $ m_{H/A}$ in the range 150--300 GeV. In order to avoid Higgs boson
conversion $A \to ZH$, Higgs boson masses are chosen with $m_A-m_H < m_Z$. It
is shown that with a proper kinematic correction applied on final state $b$-jet
four momenta, true combinations of $b$-jets can be found for simultaneous
reconstruction of both Higgs bosons through $b\bar{b}$ invariant mass
calculation. Results show that observable signals can be achieved with
statistical significance exceeding $5\sigma$ well before the target integrated
luminosity of 8 $ab^{-1}$.
