PurposeThis study aims to treat a novel system of Volterra integro-differential equations with multiple delays and variable bounds, constituting a generic numerical method based on the matrix equation and a combinatoric-parametric Charlier polynomials. The proposed method utilizes these polynomials for the matrix relations at the collocation points.Design/methodology/approachThanks to the combinatorial eligibility of the method, the functional terms can be transformed into the generic matrix relations with low dimensions, and their resulting matrix equation. The obtained solutions are tested with regard to the parametric behaviour of the polynomials with $\alpha$, taking into account the condition number of an outcome matrix of the method. Residual error estimation improves those solutions without using any external method. A calculation of the residual error bound is also fulfilled.FindingsAll computations are carried out by a special programming module. The accuracy and productivity of the method are scrutinized via numerical and graphical results. Based on the discussions, one can point out that the method is very proper to solve a system in question.Originality/valueThis paper introduces a generic computational numerical method containing the matrix expansions of the combinatoric Charlier polynomials, in order to treat the system of Volterra integro-differential equations with multiple delays and variable bounds. Thus, the method enables to evaluate stiff differential and integral parts of the system in question. That is, these parts generates two novel components in terms of unknown terms with both differentiated and delay arguments. A rigorous error analysis is deployed via the residual function. Four benchmark problems are solved and interpreted. Their graphical and numerical results validate accuracy and efficiency of the proposed method. In fact, a generic method is, thereby, provided into the literature.