This work studies the self-interacting dark matter (SIDM) scenario in the general NMSSM and beyond, where the dark matter is a Majorana fermion and the force mediator is a scalar boson. An improved analytical estimation for the dark matter (DM) self-interacting cross section is proposed, which matches well with the numerical result within O(0.1) for a wide parameter space. Due to the large couplings and light mediator in SIDM scenario, the DM/mediator will go through multiple branchings if they are produced with high energy. Based on the Monte Carlo simulation of the showers in DM sector, we obtain the multiplicities and the spectra of the DM/mediator from the Higgsino production and decay at the LHC for our benchmark points.