question: We use the total number of individuals involved in the coronavirus disease-2019 (COVID-19), namely, N, inside a specific region as a parameter in the susceptible-infected-quarantined-recovery (SIQR) model of Odagaki. Public data on the number of newly detected individuals are fitted by the numerical results of the SIQR model with optimized parameters. As a result of the optimization, we can determine the total number of individuals involved in COVID-19 inside a specific region and call such an SIQR model with a realistic total number of people involved the SIQR-N model. We then propose two methods to simulate multiple epidemic waves (MEWs), which appear in the time evolution of the number of the newly detected individuals. One is a decomposition of MEWs into independent epidemic waves that can be approximated by multiple time-derivative logistic functions (MTLF). Once the decomposition of the MEWs is completed, we fit the solution of the SIQR-N model to each MTLF using optimized parameters. Finally, we superpose the solutions obtained by multiple SIQR-N (MSIQR-N) models with the optimized parameters to fit the MEWs. The other is a set of N in the SIQR-N model as a function of time, namely, N(t), now called the SIQR-N ( t ) model. Numerical results indicate that a logistic functional approximation of N(t) fits MEWs with good accuracy. Finally, we confirm the availability of the MSIQR-N model with effects of vaccination using the recent data in Israel.