In this study, we propose a suitable solution for a non-linear system of ordinary differential equations (ODE) of the first order with the initial value problems (IVP) that contains multi variables and multi-parameters with missing real data. To solve the mentioned system, a new modified numerical simulation method is created for the first time which is called Mean Latin Hypercube Runge-Kutta (MLHRK). This method can be obtained by combining the Runge-Kutta (RK) method with the statistical simulation procedure which is the Latin Hypercube Sampling (LHS) method. The present work is applied to the influenza epidemic model in Australia in 1919  for a previous study. The comparison between the numerical and numerical simulation results is done, discussed and tabulated. The behavior of subpopulations is shown graphically. MLHRK method can reduce the number of numerical iterations of RK, and the number of LHS simulations, thus it saves time, effort, and cost.  As well as it is a faster simulation over the distribution of the LHS. The MLHRK method has been proven to be effective, reliable, and  convergent to solve a wide range of linear and nonlinear problems. The proposed method can predict the future behavior of the population under study in analyzing the behavior of some epidemiological models.