  We calculate the double charmonium production cross section within the
framework of $4\times 4$ Bethe-Salpeter Equation in the electron-positron
annihilation, at center of mass energy $\sqrt{s}= 10.6$GeV, that proceeds
through the exchange of a single virtual photon. In this calculation, we make
use of the full Dirac structure of 4D BS wave functions of these charmonia,
with the incorporation of all the Dirac covariants (both leading and
sub-leading). The calculated cross-sections for the double charmonium
productions for final states, ($J/\Psi, \eta_c$); ($\Psi',\eta_c$) ; ($J/\Psi,
\eta_c'$); and ($\Psi',\eta_c'$) are close to experimental data, and in broad
agreement with results of other theoretical models.
