We consider the propagation of surface water waves described by the Boussinesq system. Following (Molinet et al. in Nonlinearity 34:744–775, 2021), we introduce a regularized Boussinesq system obtained by adding a non-local pseudo-differential operator define by $\widehat{g_{\lambda }[\zeta ]}=|k|^{\lambda }\hat{\zeta }_{k}$ with $\lambda \in ]0,2]$ . In this paper, we display a twofold approach: first, we study theoretically the existence of an asymptotic expansion for the solution to the Cauchy problem associated to this regularized Boussinesq system with respect to the regularizing parameter $\epsilon $ . Then, we compute numerically the function coefficients of the expansion (in $\epsilon $ ) and verify numerically the validity of this expansion up to order 2. We also check the numerical $L^{2}$ stability of the numerical algorithm.