EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 1, Article Number 5843
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

An Efficient Adaptive Time-stepping Method for the
Modeling of Epidemic Dynamics

Junseok Kim

Department of Mathematics, Korea University, Seoul 02841, Republic of Korea

Abstract. In this study, we present a novel and efficient adaptive time-stepping method for mod-
eling epidemic dynamics. Examples of mathematical epidemic models include the susceptible-
infected-recovered (SIR) model, the susceptible-exposed-infected-recovered (SEIR) model, the sus-
ceptible -infected-susceptible (SIS) model, the susceptible-infected-recovered-susceptible (SIRS)
model, and the susceptible-infected-quarantined-recovered (SIQR) model. More complex models
include the maternal immunity susceptible-infected-recovered (MSIR) model, the age-structured
SEIR model, and stochastic epidemic models. These models are designed to capture specific dis-
ease characteristics, such as latency, immunity duration, and intervention impacts, and are essential
tools for studying the dynamics of infectious diseases in diverse populations. The proposed adap-
tive time-stepping method is based on the total magnitude of the summation of compartment
population differences after a single time step. Unlike other adaptive methodologies, the proposed
algorithm does not require recalculations to satisfy a given tolerance and achieves the desired accu-
racy with a single update. Therefore, the adaptive time-stepping method is both straightforward
and efficient. Several numerical tests are conducted to demonstrate the superior performance of
the proposed method.

2020 Mathematics Subject Classifications: 65L05, 65L06, 92D30

Key Words and Phrases: adaptive time-stepping, epidemic mathematical model, numerical
method, SIR model

1. Introduction

It is indisputable that mathematical epidemic models hold paramount importance
in the formulation of strategies in order to prevent the transmission of infectious dis-
eases, including COVID-19, because they provide quantitative insights into transmission
patterns and inform evidence-based policymaking [14]. These mathematical epidemic
models provide essential tools for the estimation of important epidemiological param-
eters, such as the basic reproduction number, and enable researchers to assess diverse
intervention measures [2]. This approach enables public health authorities to implement
timely and effective responses. Such measures reduce the impact of large-scale outbreaks

DOI: https://doi.org/10.29020/nybg.ejpam.v18i1.5843

Email address: cfdkim@korea.ac.kr (J. Kim)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 2 of 10

[15]. Therefore, epidemic models are indispensable in epidemiology and public health
research because they expedite a systematic examination of disease transmission and con-
trol strategies. Mathematical frameworks such as the susceptible-infected-recovered (SIR)
[18], susceptible-exposed-infected-recovered (SEIR) [10], susceptible-infected-susceptible
(SIS) [17], susceptible-infected-recovered-susceptible (SIRS) [11], susceptible-infected -
quarantined-recovered (SIQR) [13], maternal immunity susceptible-infected -recovered
(MSIR) [16], age-structured SEIR [4], and stochastic variants [12] frequently appear in
studies that resolve diverse aspects of infectious diseases. Each framework isolates essen-
tial characteristics, such as latency periods, immunity duration, or the impact of interven-
tions. These characteristics are crucial for the design of effective prevention and mitigation
strategies.

A wide range of computational methods exists for the simulation of epidemic models.
Among these approaches, robust numerical methods offer accurate and reliable predictions.
However, many standard methods rely on fixed time steps that do not adapt to significant
changes in a system’s dynamics. This shortcoming can produce unnecessary computational
overhead or diminished accuracy when the model enters regions of greater complexity or
rapid change. Consequently, there is a growing need for computational methods that can
automatically adjust time steps while maintaining efficiency.

Adaptive time-stepping approaches resolve this need by modifying the step size in
response to variations in the underlying model equations [9]. This flexibility reduces re-
dundant calculations and preserves the accuracy of the results. However, many existing
adaptive schemes require repeated recalculations to enforce a prescribed tolerance, which
diminishes their overall efficiency. This paper introduces a novel adaptive time-step se-
lection method for epidemic models that uses the total magnitude of the summation of
compartment population differences after a single step as the sole criterion for time-step
adjustment. This approach achieves the required accuracy with a single update and elim-
inates additional recalculations to satisfy tolerance requirements. Furthermore, the crite-
rion for adaptive time-stepping is more intuitive compared to previous criteria based on
error estimation. The proposed method is straightforward to implement and compatible
with the aforementioned broad range of epidemic models. Numerical experiments confirm
the method’s efficiency.

The organization of this article is described as follows. Section 2 presents the proposed
adaptive time-stepping method. Section 3 provides computational tests to validate the
efficiency and accuracy of the proposed method. Section 4 concludes the article with a
discussion of the main findings and potential future work. Additionally, the Appendix
includes the MATLAB source code for the proposed method, which can be used as a
computational tool for further analysis and replication of results.

2. Proposed Computational Method

The standard SIR model is considered to demonstrate the proposed adaptive time-
stepping algorithm for mathematical epidemic models [1]. The SIR model describes the
transmission of infectious diseases by dividing the population into three groups: susceptible



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 3 of 10

(S), representing individuals at risk of infection; infected (I), representing those currently
infected and capable of transmitting the disease; and recovered (R), representing individ-
uals who have recovered and gained immunity. The standard SIR model is formulated as
follows:

dS(t)

dt
= −βS(t)I(t), (1)

dI(t)

dt
= βS(t)I(t)− γI(t), (2)

dR(t)

dt
= γI(t). (3)

Here, S(t), I(t), and R(t) represent the numbers of susceptible, infected, and recovered
individuals at time t, corresponding to those not currently infected, those actively infected,
and those who have recovered, respectively. The parameter β is the average number of
random contacts an individual has per unit time, while 1/γ denotes the average duration
for an infected individual to recover [3]. To effectively solve the governing equations of
the SIR model (1)–(3) and capture the dynamics of disease transmission, we propose a
novel numerical method that uses an adaptive time-stepping approach designed to the
nonuniform temporal changes in the model. Let Sn = S(tn), In = I(tn), and Rn = R(tn)
for n = 1, . . . , where tn = tn−1 + ∆tn, for n ≥ 1, ∆tn is the variable time step size,
and t0 = 0. The SIR model (1)–(3) can be discretized using a finite difference method [8]
as follows:

Sn+1 − Sn

∆tn+1
= −βSnIn, (4)

In+1 − In
∆tn+1

= βSnIn − γIn, (5)

Rn+1 −Rn

∆tn+1
= γIn, (6)

which can be rewritten as follows:

Sn+1 = Sn −∆tn+1βSnIn, (7)

In+1 = In +∆tn+1(βSnIn − γIn), (8)

Rn+1 = Rn +∆tn+1γIn. (9)

Using a recently developed methodology [6], we describe an adaptive time-stepping al-
gorithm for epidemic mathematical models. We define a tolerance tol to set the maximum
allowable displacement of the computational solution during a time step. By multiplying
both sides of Eqs. (4)–(6) by ∆tn+1, the following equations are derived.

Sn+1 − Sn = −∆tn+1βSnIn, (10)

In+1 − In = ∆tn+1(βSnIn − γIn), (11)

Rn+1 −Rn = ∆tn+1γIn. (12)



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 4 of 10

Let evol be the total magnitude of the summation of each compartment’s population
differences after a single time step:

evol = |Sn+1 − Sn|+ |In+1 − In|+ |Rn+1 −Rn| (13)

= ∆tn+1βSnIn + |∆tn+1(βSnIn − γIn)|+∆tn+1γIn

= ∆tn+1(βSnIn + |βSnIn − γIn|+ γIn).

In this study, it is important to note that we use Eq. (13) for evol, which indicates that
our objective is to control the evolution by accounting for population differences among
the considered groups. We may use other criterion such as

evol =
√

(Sn+1 − Sn)2 + (In+1 − In)2 + (Rn+1 −Rn)2. (14)

Let us define tol as a given tolerance. We require that ∆tn+1 satisfy the following
condition: evol ≤ tol, that is,

∆tn+1(βSnIn + |βSnIn − γIn|+ γIn) ≤ tol. (15)

From Eq. (15), we obtain the following constraint on the time step sizes:

∆tn+1 ≤ tol

βSnIn + |βSnIn − γIn|+ γIn
. (16)

This constraint ensures that the time step is adjusted appropriately to satisfy the
tolerance requirement tol. Then, we use

∆tn+1 = min

(
tol

βSnIn + |βSnIn − γIn|+ γIn
, ∆tmax

)
. (17)

Here, the maximum time step size, ∆tmax, is predetermined to guarantee stability. Hence,
the numerical solutions at each time step are provided in Eqs. (10)–(12). This process is
repeated while satisfying the condition tn+∆tn+1 ≤ T to compute the numerical solutions
of S(T ), I(T ), and R(T ). When selecting adaptive time steps, it is possible for tn+1 to
exceed the specified final time T . To prevent this, if tn+1 = tn+∆tn+1 > T , the time step
tn+1 is adjusted to T − tn.

3. Computational Experiments

This section presents a series of computational experiments designed to thoroughly
evaluate the effectiveness, accuracy, and high performance of the proposed method across
different scenarios, which highlights its robustness and practical applicability in solving
epidemic problems with improved computational efficiency and reliability.



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 5 of 10

3.1. Effect of the tolerance, tol

Figure 1(a) shows the temporal evolution of the susceptible S(t), infected I(t), and
recovered R(t) populations, along with the adaptive time-stepping behavior in the discrete
time grid. The parameters are β = 0.3, γ = 0.1, and initial conditions S0 = 0.99, I0 =
0.01, and R0 = 0. S(t) decreases over time and represents the decline of the susceptible
population due to infection. I(t) initially increases, reaching a peak, and then declines
as individuals recover. R(t) increases steadily, showing the accumulation of recovered
individuals. Figure 1(b) shows the adaptive time step size (∆tn+1) and the corresponding
time grid (tn+1). The time step size adjusts dynamically, decreasing during rapid changes
(e.g., the peak of I(t)) and increasing as the system approaches a steady state. The
adaptive method ensures both accuracy and computational efficiency. In this computation,
a tolerance of tol = 0.1 and a maximum step size of ∆tmax = 3 are used.

(a)

(b)
Figure 1: (a) Computational results for S(t), I(t), and R(t) with an adaptive discrete time grid (red circles).
(b) The adaptive discrete time grid and corresponding time steps. The parameters used are β = 0.3, γ = 0.1,
S0 = 0.99, I0 = 0.01, R0 = 0, T = 120, tol = 0.1, and ∆tmax = 3.

When the tolerance is reduced to tol = 0.01, the adaptive time step size ∆tn+1 de-
creases further during rapid transitions in the system (e.g., around the peak of I(t)) to
maintain higher accuracy as shown in Fig. 2(a). As a result, more time steps are re-
quired to accurately capture the dynamics of the system. As the system approaches a
steady state, the step size increases, which optimizes computational efficiency and main-
tains accuracy. The adaptive method ensures reliable numerical results for all populations.
However, compared to the case with tol = 0.1, the computation demands higher precision
and involves more steps to satisfy the stricter tolerance requirement as shown in Fig. 2(b).



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 6 of 10

For this computation, the maximum step size of ∆tmax = 3 is used.

(a)

(b)
Figure 2: (a) Computational results for S(t), I(t), and R(t) with an adaptive discrete time grid (red circles).
(b) The adaptive discrete time grid and corresponding time steps. The parameters used are β = 0.3, γ = 0.1,
S0 = 0.99, I0 = 0.01, R0 = 0, T = 120, tol = 0.01, and ∆tmax = 3.

Figure 3 illustrates the temporal evolution of the susceptible S(t), infected I(t), and
recovered R(t) populations, alongside the adaptive time-stepping behavior for a discrete
time grid. The parameters used are β = 0.3, γ = 0.1, S0 = 0.99, I0 = 0.01, R0 = 0,
with a total simulation time of T = 120, a tolerance of tol = 0.001, and a maximum step
size of ∆tmax = 3. Figure 3(a) shows the evolution of the populations. S(t) decreases
over time due to infection, I(t) initially rises to a peak before declining as individuals
recover, and R(t) steadily increases, representing the cumulative number of recovered
individuals. Figure 3(b) presents the adaptive time step size ∆tn+1 and the corresponding
time points tn+1. Compared to results with larger tolerances, the smaller tol = 0.001
results in significantly smaller time steps during periods of rapid change, such as the peak
of I(t). As the system approaches a steady state, the step size increases and ensures
computational efficiency while maintaining accuracy. This stricter tolerance results in
finer resolution and more accurate modeling of the system dynamics, especially during
critical transitions. In this test, using a fully explicit Euler method results in a CPU time
that is double that of the adaptive algorithm.

The value of tol is selected based on the specific needs and objectives of the user and
is carefully determined according to the scale and characteristics of the problem under
consideration. This ensures that the chosen tolerance value aligns with the desired level
of accuracy and the nature of the computational model being analyzed.



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 7 of 10

(a)

(b)
Figure 3: (a) Computational results for S(t), I(t), and R(t) with an adaptive discrete time grid (red circles).
(b) The adaptive discrete time grid and corresponding time steps. The parameters used are β = 0.3, γ = 0.1,
S0 = 0.99, I0 = 0.01, R0 = 0, T = 120, tol = 0.001, and ∆tmax = 3.

3.2. Effect of the maximum step size, ∆tmax

Next, we investigate the effect of the maximum step size, ∆tmax. The parameters are
set as β = 0.3, γ = 0.1, tol = 0.1, with initial conditions S0 = 0.99, I0 = 0.01, and R0 = 0.
The computations are conducted for three different values of the maximum step size of
∆tmax. Figure 4(a) presents the computational results for S(t), I(t), and R(t) under three
different values of the maximum step size, ∆tmax. Figure 4(b) illustrates the adaptive
discrete time steps corresponding to these three values of the maximum step size. From
these results, it can be observed that the numerical solution becomes less accurate when an
excessively large ∆tmax is used. This is because, during the initial stages when the solution
is smooth, larger time steps are used, which leads to significant errors. Consequently, using
an excessively large maximum time step is not advisable.

4. Conclusions

In this study, we proposed an efficient adaptive time-stepping method for modeling epi-
demic dynamics. By dynamically adjusting the time step size based on the magnitude of
changes in population compartments, the method achieved both high accuracy and compu-
tational efficiency. Numerical experiments demonstrated the effectiveness of the proposed
approach in capturing the dynamics of susceptible, infected, and recovered populations
under various tolerance levels. The results confirmed that the adaptive method provides



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 8 of 10

(a)

(b)
Figure 4: (a) Computational results for S(t), I(t), and R(t) for three different values of the maximum step
size of ∆tmax. (b) The adaptive discrete time steps. The parameters used are β = 0.3, γ = 0.1, S0 = 0.99,
I0 = 0.01, R0 = 0, T = 120, and tol = 0.1.

finer resolution during rapid transitions, such as near the infection peak, and optimizes
step sizes as the system approaches equilibrium. The proposed method is straightforward
to implement and eliminates the need for repeated recalculations to satisfy tolerance con-
straints, making it highly efficient. These advantages make it a promising tool for a wide
range of epidemic models, including those with additional compartments or complex dy-
namics. A limitation of the proposed method lies in its reliance on a first-order scheme.
Developing a second-order scheme would be beneficial to improve accuracy, although it
poses significant challenges. For future research, the proposed method could be extended
to solve partial differential equations, such as the diffusion equation [7] and fractional
models [5], to broaden its applicability and demonstrate its effectiveness in solving more
complex problems.

Acknowledgements

This work was supported by the National Research Foundation(NRF), Korea, under
project BK21 FOUR. I sincerely thank the reviewers for their valuable comments and



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 9 of 10

constructive suggestions, which have significantly contributed to improving the quality
and clarity of this manuscript.

References

[1] I. Cooper, A. Mondal, and C.G. Antonopoulos. A SIR model assumption for
the spread of COVID-19 in different communities. Chaos, Solitons and Fractals,
139:110057, 2020.

[2] A. Gómez-Corral, F. Palacios-Rodŕıguez, and M.T. Rodŕıguez-Bernal. On the exact
reproduction number in SIS epidemic models with vertical transmission. Computa-
tional and Applied Mathematics, 42:291, 2023.

[3] Y. Hwang, S. Kwak, and J. Kim. Long-Time Analysis of a Time-Dependent SUC
Epidemic Model for the COVID-19 Pandemic. Journal of Healthcare Engineering,
2021:5877217, 2021.

[4] F.Z. Jaouimaa, D. Dempsey, S. Van Osch, S. Kinsella, K. Burke, J. Wyse, and
J. Sweeney. An age-structured SEIR model for COVID-19 incidence in Dublin, Ireland
with framework for evaluating health intervention cost. PLOS ONE, 16(12):e0260632,
2021.

[5] M. Khader. Using Modified Fractional Euler Formula for Solving the Fractional
Smoking Model. European Journal of Pure and Applied Mathematics, 17(4):2676–
2691, 2024.

[6] H. Kim, S. Kwak, M. Mohammed, S. Kang, S. Ham, and J. Kim. An efficient and ac-
curate adaptive time-stepping method for the Landau–Lifshitz equation. Algorithms,
18(1):1, 2025.

[7] J. Kim, S. Kwak, H. G. Lee, Y. Hwang, and S. Ham. A maximum principle of
the Fourier spectral method for diffusion equations. Electronic Research Archive,
31(9):5396–5405, 2023.

[8] S. Kwak, S. Kang, S. Ham, Y. Hwang, G. Lee, and J. Kim. An unconditionally stable
difference scheme for the two-dimensional modified Fisher–Kolmogorov–Petrovsky–
Piscounov equation. Journal of Mathematics, 2023:Article ID 5527728, 2023.

[9] C. Lee, J. Park, S. Kwak, S. Kim, Y. Choi, S. Ham, and J. Kim. An adaptive
time-stepping algorithm for the Allen–Cahn equation. Journal of Function Spaces,
2022:Article ID 2731593, 2022.

[10] Q. Liu and D. Jiang. Stationary distribution and probability density for a stochastic
SEIR-type model of coronavirus (COVID-19) with asymptomatic carriers. Chaos,
Solitons and Fractals, 169:113256, 2023.

[11] F. Nill. Endemic oscillations for SARS-COV-2 Omicron—A SIRS model analysis.
Chaos, Solitons and Fractals, 173:113678, 2023.

[12] M. Pájaro, N.M. Fajar, A.A. Alonso, and I. Otero-Muras. Stochastic SIR model
predicts the evolution of COVID-19 epidemics from public health and wastewater
data in small and medium-sized municipalities: A one year study. Chaos, Solitons
and Fractals, 164:112671, 2022.

[13] C.M. Pinto, J.A. Tenreiro Machado, and C. Burgos-Simón. Modified SIQR model for



J. Kim / Eur. J. Pure Appl. Math, 18 (1) (2025), 5843 10 of 10

the COVID-19 outbreak in several countries. Mathematical Methods in the Applied
Sciences, 47(5):3273–3288, 2024.

[14] A.K. Sabherwal, A. Sood, and M.A. Shah. Evaluating mathematical models for
predicting the transmission of COVID-19 and its variants towards sustainable health
and well-being. Discover Sustainability, 5:38, 2024.

[15] S. Sarin, S. Singh, S. Kumar, S. Goyal, B. Gupta, V. Arya, and K. Chui. SEIR-
driven semantic integration framework: Internet of Things-enhanced epidemiological
surveillance in COVID-19 outbreaks using recurrent neural networks. IET Cyber-
Physical Systems: Theory and Applications, 9:135–149, 2024.

[16] F.J. Schwarzendahl, J. Grauer, B. Liebchen, and H. Löwen. Mutation induced infec-
tion waves in diseases like COVID-19. Scientific Reports, 12(1):9641, 2022.

[17] A.A. Thirthar, H. Abboubakar, A. Khan, and T. Abdeljawad. Mathematical modeling
of the COVID-19 epidemic with fear impact. AIMS Mathematics, 8(3):6447–6465,
2023.

[18] Y. Wang, F. Yuan, Y. Song, H. Rao, L. Xiao, H. Guo, et al. Prediction of cross-border
spread of the COVID-19 pandemic: A predictive model for imported cases outside
China. PLOS ONE, 19(4):e0301420, 2024.