In this article, a numerical method integrated with statistical data simulation technique is introduced to solve a nonlinear system of ordinary differential equations with multiple random variable coefficients. The utilization of Monte Carlo simulation with central divided difference formula of finite difference (FD) method is repeated n times to simulate values of the variable coefficients as random sampling instead being limited as real values with respect to time. The mean of the n final solutions via this integrated technique, named in short as mean Monte Carlo finite difference (MMCFD) method, represents the final solution of the system. This method is proposed for the first time to calculate the numerical solution obtained for each subpopulation as a vector distribution. The numerical outputs are tabulated, graphed, and compared with previous statistical estimations for 2013, 2015, and 2030, respectively. The solutions of FD and MMCFD are found to be in good agreement with small standard deviation of the means, and small measure of difference. The new MMCFD method is useful to predict intervals of random distributions for the numerical solutions of this epidemiology model with better approximation and agreement between existing statistical estimations and FD numerical solutions.