This paper introduces a non-conventional approach with multi-dimensional random sampling to solve a cocaine abuse model with statistical probability. The mean Latin hypercube finite difference (MLHFD) method is proposed for the first time via hybrid integration of the classical numerical finite difference (FD) formula with Latin hypercube sampling (LHS) technique to create a random distribution for the model parameters which are dependent on time [Formula: see text]. The LHS technique gives advantage to MLHFD method to produce fast variation of the parameters’ values via number of multidimensional simulations (100, 1000 and 5000). The generated Latin hypercube sample which is random or non-deterministic in nature is further integrated with the FD method to complete one cycle of LHS-FD simulation iteration. This process is repeated until [Formula: see text] final iterations of LHS-FD are obtained. The means of these [Formula: see text] final solutions (MLHFD solutions) are tabulated, graphed and analyzed. The numerical simulation results of MLHFD for the SEIR model are presented side-by-side with deterministic solutions obtained from the classical FD scheme and homotopy analysis method with Pade approximation (HAM-Pade). The present MLHFD results are also compared with the previous non-deterministic statistical estimations from 1995 to 2015. Good agreement between the two is perceived with small errors. MLHFD method can be used to predict future behavior, range and prediction interval for the epidemic model solutions. The expected profiles of the cocaine abuse subpopulations are projected until the year 2045. Both the statistical estimations and the deterministic results of FD and HAM-Pade are found to be within the MLHFD prediction intervals for all the years and for all the subpopulations considered.