EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6007 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Numerical Solution for Nash Differential Games Based on the Runge Kutta 4th-Order Method Abd El-Monem A. Megahed1,∗, Nesreen M. Kamel1, I. M. Hanafy2, Nora A. Omar2 1 Department of Basic Science, Faculty of Computers and Informatics, Suez Canal University, Ismailia 41522, Egypt 2 Department of Mathematics and Computer Science, Faculty of Science, Port Said University, Egypt Abstract. In this paper, we present a numerical solution for an open-loop Nash differential game modeling competition between two firms. Using the fourth-order Runge-Kutta method, we com- puted the numerical solution and analyzed the stability of the open-loop Nash equilibrium. Addi- tionally, we examined the uniform convergence of the solution. Finally, an illustrative example is presented to clarify the results, accompanied by figures to illustrate the findings. 2020 Mathematics Subject Classifications: 91A23, 49N05, 49N70, 49N90 Key Words and Phrases: Open-loop Nash differential game, Runge-Kutta fourth-order method, stability, uniform convergence 1. Introduction Differential games are a branch of game theory that focus on studying and develop- ing optimal control strategies for dynamic systems influenced by the decisions of multiple players [1]. An open-loop Nash differential game is a specific type of differential game in which the strategies of the players depend solely on time (t) rather than the current state of the dynamic system. In other words, players predefine their strategies at the start of the game and do not update them based on subsequent observations or system states. In game theory, differential games are used to model and analyze conflicts, such as competition, within dynamical systems. Differential equations play a key role in many applications in physics, engineering, and the modeling of natural phenomena [2]. An effective approach for solving differential games is through numerical methods. For ex- ample, Dehghan Banadaki and Navidi [3] studied open-loop Nash differential games using ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6007 Email addresses: abdelmoneam mgahed@ci.suez.edu.eg (A. A. Megahed), nesreenmostafakamelbs@ci.suez.edu.eg (N. M. Kamel), ihanafy@hotmail.com (I. M. Hanafy), noraalaa19862013@gmail.com (N. A. Omar) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 2 of 22 the Legendre Tau method combined with the fourth-order Runge–Kutta method. Ku- vshinov and Osipov [4] investigated Stackelberg solutions for linear positional differential games by applying the polyhedron method. Megahed et al. [5] studied an open-loop Nash differential game and employed the Picard method to approximate the solution of a model describing competition between two firms, incorporating market share, advertising efforts, and advertising costs. Megahed et al. [6] studied a zero-sum differential game for model- ing coronavirus dynamics, solving it using the homotopy perturbation method combined with a new iterative approach. The Stackelberg differential game of E-differentiable and E-convex functions was applied to combat terrorism, considering government actions [7]. Hemeda [8] introduced the Integral Iterative Method (IIM), a modification of the Picard method, as a numerical technique for solving nonlinear integro-differential equations and systems. Megahed et al. [9] investigated an open-loop Nash differential game, where they derived the necessary conditions for equilibrium and analyzed the existence and uniform convergence of the solutions using the Picard method. Illustrative figures were provided to demonstrate applications in economic, financial, and industrial contexts. Kassem et al. [10] discussed a Nash-collaborative approach for a differential game involv- ing multiple governments and terrorist organizations, analyzing government cooperation and the role of each government in counterterrorism. Youness [11] introduced a differ- ential game termed “Nash coalitional,” which extends the traditional Nash equilibrium framework by incorporating partial cooperation among players, enabling them to achieve mutually beneficial outcomes while still accounting for individual objectives. Youness et al. [12] examined the necessary conditions for determining optimal strategies in fuzzy continuous differential games under Nash equilibrium. Sun et al. [13] introduced a lin- ear–quadratic stochastic two-person nonzero-sum differential game with both open-loop and closed-loop Nash equilibria. Engwerda [14] analyzed the open-loop Nash equilib- rium in a linear-quadratic (LQ) differential game. Min-max differential games with fuzzy objectives and controls, as well as large-scale differential games, have been discussed in [15, 16]. Megahed [17] analyzed terrorism dynamics through a two-player differential game involving the International Terrorism Organization (ITO). Youness et al. [18] examined a min-max differential game whose state trajectory is governed by a Cauchy initial value problem (CIVP). They provided both analytical and approximate numerical solutions us- ing the Picard method, demonstrating the model’s effectiveness in dynamic optimization contexts. In this paper, we present the numerical solution of an open-loop Nash differential game using the fourth-order Runge–Kutta method [19]. The system dynamics, representing competition between two firms in the market, are modeled by differential equations. Each player aims to optimize their objective while accounting for the impact of both their own actions and those of the other player. Section 2 formulates the dynamical system of the problem and derives the necessary conditions for an open-loop Nash equilibrium. Section 3 presents the numerical solution using the fourth-order Runge–Kutta method and discusses its convergence. Section 4 analyzes the stability of the numerical solution, and Section 5 concludes the paper. A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 3 of 22 2. Problem Formulation In this section, we discuss the solution of an open-loop Nash differential game between two firms over the time interval t ∈ [0, T ]. The system dynamics can be expressed as follows: dx dt = u1(t) ( 1− x(t) ) − u2(t)x(t), (1) where t ∈ [0, T ], x(0) = x0, subject to the constraint 0 ≤ x(t) ≤ 1. (2) where x(t) denotes the market share of Firm 1 at time t, and 1 − x(t) represents the market share of Firm 2. The control variables u1(t) and u2(t) are defined as follows: u1(t) corresponds to the advertising efforts of Firm 1 at time t, while u2(t) corresponds to the advertising efforts of Firm 2 at time t. For Firm 1, the number of customers increases due to its own advertising efforts, whereas the advertising efforts of Firm 2 decrease the number of customers of Firm 1. The state dynamics and the payoff functionals are described as follows: Ji ( u1(t), u2(t) ) = ∫ T 0 Ii ( x(t), u1(t), u2(t), t ) dt, i = 1, 2. (3) The payoff functionals of the two firms in problems (1) and (2) are defined as follows: J1(u1(t), u2(t)) = ∫ T 0 e−r1t[ϕ1x(t)− C1u1(t)] dt, J2(u1(t), u2(t)) = ∫ T 0 e−r2t[ϕ2(1− x(t))− C2u2(t)] dt. (4) where ri is the interest rate of Firm i, ϕi is the fractional revenue potential of Firm i, and Ci(s) denotes the advertising cost function. We assume that Ci(s) = Pi 2 s2, Pi > 0, i = 1, 2. f(x(t), u1(t), u2(t), t) = u1(t) ( 1− x(t) ) − u2(t)x(t), I1(x(t), u1(t), u2(t), t) = ϕ1 x(t)− C1 ( u1(t) ) = ϕ1 x(t)− P1 2 u1(t) 2, I2(x(t), u1(t), u2(t), t) = ϕ2 (1− x(t))− C2 ( u2(t) ) = ϕ2 (1− x(t))− P2 2 u2(t) 2. (5) Theorem 1. Let f(x(t), u1(t), u2(t), t) and Ii(x(t), u1(t), u2(t), t) be continuously differ- entiable functions on Rn, i = 1, 2. Specifically, f(x(t), u1(t), u2(t), t) : Rn × Rs × [0, T ] → R, f ∈ C1, A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 4 of 22 where s = ∑N j=1 sj, i ̸= j, and Ii(x(t), u1(t), u2(t), t) : Rn × Rs × [0, T ] → R, Ii ∈ C1, i = 1, 2. If u∗i (t), 0 ≤ t ≤ T , is an open-loop Nash equilibrium solution, and x∗(t), 0 ≤ t ≤ T , is the corresponding state trajectory. Then, there exist two costate vectors λi(t) : [0, T ] → Rn and two Hamiltonian functions . Hi ( λi(t), x(t), u1(t), u2(t), t ) = Ii ( x(t), u1(t), u2(t), t ) +λi(t) ⊤f ( x(t), u1(t), u2(t), t ) , i = 1, 2. (6) such that the following necessary conditions hold: dx∗(t) dt = f ( x∗(t), u∗1(t), u ∗ 2(t), t ) , x∗(0) = x0, dλ∗ 1(t) dt = − ∂H1 ( λ1(t), x ∗(t), u∗1(t), u ∗ 2(t), t ) ∂x , dλ∗ 2(t) dt = − ∂H2 ( λ2(t), x ∗(t), u∗1(t), u ∗ 2(t), t ) ∂x , ∂H1 ( λ1(t), x ∗(t), u∗1(t), u ∗ 2(t), t ) ∂u1 = 0, ∂H2 ( λ2(t), x ∗(t), u∗1(t), u ∗ 2(t), t ) ∂u2 = 0. (7) with initial and terminal conditions x∗(0) = x0, λ1(T ) = 0, λ2(T ) = 0. (8) The proof of this theorem is presented in [20]. 3. The Numerical Solution by Using the Fourth-Order Runge-Kutta Method 3.1. Existence and Convergence of the Solution Existence of the Solution After applying the necessary conditions for an open-loop Nash equilibrium differential game, the problem (1)–(5) is reduced to the following Hamiltonian functions. For Player 1, the Hamiltonian is defined as H1(λ1(t), x(t), u1(t), u2(t), t) = ϕ1x(t)− P1 2 u21(t)+λ1(t) T [ u1(t)(1−x(t))−u2(t)x(t) ] . (9) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 5 of 22 And for Player 2, the Hamiltonian is H2(λ2(t), x(t), u1(t), u2(t), t) = ϕ2 ( 1−x(t) ) − P2 2 u22(t)+λ2(t) T [ u1(t)(1−x(t))−u2(t)x(t) ] . (10) dx dt = u1(t) ( 1− x(t) ) − u2(t)x(t), x(0) = x0, dλ1(t) dt = λ1(t) [ u1(t) + u2(t) ] − ϕ1, λ1(T ) = 0, dλ2(t) dt = λ2(t) [ u1(t) + u2(t) ] + ϕ2, λ2(T ) = 0, u1(t) = λ1(t) P1 ( 1− x(t) ) , u2(t) = −λ2(t) P2 x(t), (ci)n = pi 2 ( u2i ) n , i = 1, 2. (11) The existence of a solution to system (11) is discussed in [9]. Uniform Convergence of the Solution Now, we apply the fourth-order Runge–Kutta method to discuss the uniform conver- gence of the solution. The scheme is given by xn+1 = xn + h 6 Φ1(xn, λ1,n, λ2,n, tn), λ1,n+1 = λ1,n + h 6 Φ2(xn, λ1,n, λ2,n, ϕ1, tn), λ2,n+1 = λ2,n + h 6 Φ3(xn, λ1,n, λ2,n, ϕ2, tn), (12) where Φ1(xn, λ1,n, λ2,n, tn) = K1x + 2K2x + 2K3x +K4x, Φ2(xn, λ1,n, λ2,n, ϕ1, tn) = K1λ1 + 2K2λ1 + 2K3λ1 +K4λ1 , Φ3(xn, λ1,n, λ2,n, ϕ2, tn) = K1λ2 + 2K2λ2 + 2K3λ2 +K4λ2 . (13) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 6 of 22 The Runge-Kutta variables ki are defined as K1x = hFx(xn, tn), K2x = hFx ( xn + h 2 , tn + K1x 2 ) , K3x = hFx ( xn + h 2 , tn + K2x 2 ) , K4x = hFx (xn + h, tn +K3x) , K1λ1 = hFλ1(xn, tn), K2λ1 = hFλ1 ( xn + h 2 , tn + K1λ1 2 ) , K3λ1 = hFλ1 ( xn + h 2 , tn + K2λ1 2 ) , K4λ1 = hFλ1 (xn + h, tn +K3λ1) , K1λ2 = hFλ2(xn, tn), K2λ2 = hFλ2 ( xn + h 2 , tn + K1λ2 2 ) , K3λ2 = hFλ2 ( xn + h 2 , tn + K2λ2 2 ) , and K4λ2 = hFλ2 (xn + h, tn +K3λ2) . (14) The sequence xn+1, λ1,n+1, and λ2,n+1 can be written as the finite series xn+1 = x0 + n∑ j=0 (xj+1 − xj) , λ1,n+1 = λ1,0 + n∑ j=0 (λ1,j+1 − λ1,j) , λ2,n+1 = λ2,0 + n∑ j=0 (λ2,j+1 − λ2,j) . (15) If xn+1, λ1,n+1 and λ2,n+1 are convergent, then the infinite series ∞∑ j=0 (xj+1 − xj), ∞∑ j=0 ( λ1,j+1 − λ1,j ) , ∞∑ j=0 (( λ2,j+1 − λ2,j )) are convergent, and the solution will be x, λ1 and λ2. lim n→∞ xn+1 = x, lim n→∞ λ1,n+1 = λ1, lim n→∞ λ2,n+1 = λ2. (16) Assume that 1. The functions Φi : Rn × Rn × Rn × [0, T ] → R, i = 1, 2, 3, A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 7 of 22 are continuous, and there exist positive constants Mi such that |Φi| ≤ Mi, i = 1, 2, 3. 2. Each function Φi satisfies the Lipschitz condition with Lipschitz constants Li, where 0 < Li < 1, i = 1, 2, 3, such that |Φi(z1)− Φi(z2)| ≤ Li∥z1 − z2∥, ∀z1, z2 ∈ Rn × Rn × Rn × [0, T ]. In particular, |Φ1(x, λ1, λ2, t)− Φ1(y, λ1, λ2, t)| < L1|x− y|, (17) |Φ2(x, λ1, λ2, t)− Φ2(x, α, λ2, t)| < L2|λ1 − α|, (18) |Φ3(x, λ1, λ2, t)− Φ3(x, λ1, γ, t)| < L3|λ2 − γ|. (19) If the three series converge, then the three sequences xn+1, λ1,n+1 and λ2,n+1 will converge to x, λ1 and λ2, respectively. To discuss the uniform convergence of xn+1, λ1,n+1 and λ2,n+1, we consider the following three associated series: ∞∑ n=0 (xn+1 − xn), ∞∑ n=0 (λ1,n+1 − λ1,n), ∞∑ n=0 (λ2,n+1 − λ2,n). (20) For n = 0, x1 − x0 = h 6 Φ1(xn, λ1,n, λ2,n, t), |x1 − x0| = ∣∣∣∣h6Φ1(xn, λ1,n, λ2,n, t) ∣∣∣∣ , |x1 − x0| ≤ h 6 M1. (21) Similarly, |λ1,1 − λ1,0| ≤ h 6 M2, |λ2,1 − λ2,0| ≤ h 6 M3. (22) Now, we will get an estimation for(xn+1 − xn), (λ1,n+1 − λ1,n) and (λ2,n+1 − λ2,n) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 8 of 22 |xn+1 − xn| = ∣∣∣∣(xn − xn−1) + h 6 (Φ1,n − Φ1,n−1) ∣∣∣∣ , |xn+1 − xn| ≤ |xn − xn−1|+ h 6 |Φ1,n − Φ1,n−1| , |xn+1 − xn| ≤ |xn − xn−1|+ h 6 |L1(xn − xn−1)| , |xn+1 − xn| ≤ (1 + L1 h 6 ) |xn − xn−1| . (23) Similarly, |λ1,n+1 − λ1,n| ≤ (1 + h 6 L2) |λ1,n − λ1,n−1| , |λ2,n+1 − λ2,n| ≤ (1 + h 6 L3) |λ2,n − λ2,n−1| . (24) At n = 1, |x2 − x1| ≤ (1 + L1 h 6 ) |x1 − x0| |x2 − x1| ≤ h 6 M1(1 + L1 h 6 ). (25) At n = 2, |x3 − x2| ≤ (1 + L1 h 6 ) |x2 − x1| , |x3 − x2| ≤ h 6 M1(1 + L1 h 6 )2. (26) At n = 3, |x4 − x3| ≤ (1 + L1 h 6 ) |x3 − x2| , |x4 − x3| ≤ h 6 M1(1 + L1 h 6 )3. (27) And so on |xn+1 − xn| ≤ h 6 M1(1 + L1 h 6 )n. (28) Similarly, |λ1,n+1 − λ1,n| ≤ h 6 M2(1 + L1 h 6 )n, |λ2,n+1 − λ2,n| ≤ h 6 M2(1 + L1 h 6 )n. (29) Since Li ≤ 1, Mi ≤ 1 and h ≤ 1 i = 1, 2, 3. Then the series ∑∞ n=0 (xn+1 − xn) , ∑∞ n=0 (λ1,n+1 − λ1,n) , ∑∞ n=0 (λ2,n+1 − λ2,n) , and the sequences xn+1, λ1,n+1, λ2,n+1 are uniformly convergent . Since Li ≤ 1, Mi ≤ 1, and h ≤ 1 for i = 1, 2, 3, the series ∑∞ n=0(xn+1 − xn),∑∞ n=0(λ1,n+1 − λ1,n), and ∑∞ n=0(λ2,n+1 − λ2,n), as well as the sequences xn+1, λ1,n+1, and λ2,n+1, are uniformly convergent. A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 9 of 22 3.2. Numerical Solution We now apply the fourth-order Runge–Kutta method to obtain the numerical solution for the problem (11). By substituting the controls u1 and u2 into the state and costate equations, we obtain: dx dt = λ1 P1 (1− x)2 + λ2 P2 x2, x(0) = x0, dλ1(t) dt = λ2 1 P1 (1− x)− λ1λ2(t) P2 x− ϕ1, λ1(T ) = 0, dλ2(t) dt = λ1λ2(t) P1 (1− x)− λ2 2 P2 x+ ϕ2, λ2(T ) = 0. (30) Then, the system becomes : Fx(xn, tn) = (λ1(t))n P1 (1− xn(t)) 2 + (λ2(t))n P2 x2n, Fλ1(xn, tn) = (λ1(t) 2)n P1 (1− xn(t))− (λ1(t))n(λ2(t))n P2 xn(t)− ϕ1, Fλ2(xn, tn) = (λ1(t))n(λ2(t))n P1 (1− xn(t))− (λ2 2(t))n P2 xn(t) + ϕ2, (ci)n+1 = pi 2 (u2i )n+1, i = 1, 2. (31) with the initial and terminal conditions: x(0) = x0, (λ1)0 = 0, (λ2(T ))0 = 0. The general form after applying the Runge-Kutta 4th-order method is as follows, xn+1 = xn + h 6 (K1x + 2K2x + 2K3x +K4x), λ1,n+1 = λ1,n + h 6 (K1λ1 + 2K2λ1 + 2K3λ1 +K4λ1), λ2,n+1 = λ2,n + h 6 (K1λ2 + 2K2λ2 + 2K3λ2 +K4λ2). (32) Where A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 10 of 22 K1x = hFx(xn, tn), K2x = hFx ( xn + h 2 , tn + K1x 2 ) , K3x = hFx ( xn + h 2 , tn + K2x 2 ) , K4x = hFx (xn + h, tn +K3x) , K1λ1 = hFλ1(xn, tn), K2λ1 = hFλ1 ( xn + h 2 , tn + K1λ1 2 ) , K3λ1 = hFλ1 ( xn + h 2 , tn + K2λ1 2 ) , K4λ1 = hFλ1 (xn + h, tn +K3λ1) , K1λ2 = hFλ2(xn, tn), K2λ2 = hFλ2 ( xn + h 2 , tn + K1λ2 2 ) , K3λ2 = hFλ2 ( xn + h 2 , tn + K2λ2 2 ) , and K4λ2 = hFλ2 (xn + h, tn +K3λ2) . (33) At n=0, x1 = x0 + h 6 (K1x + 2K2x + 2K3x +K4x) (34) Where K1x = hFx(x0, t0) = h [ λ1,0(t) P1 (1− x0(t)) 2 + (λ2(t))0 P2 x20 ] = 0, K2x = hFx ( x0 + h 2 , t0 + K1x 2 ) = h [ λ1(t0 + K1x 2 ) P1 (1− (x0 + h 2 ))2 + λ2(t0 + K1x 2 ) P2 (x0 + h 2 )2 ] = 0 K3x = hFx ( x0 + h 2 , t0 + K2x 2 ) = h [ λ1(t0 + K2x 2 ) P1 (1− (x0 + h 2 ))2 + λ2(t0 + K2x 2 ) P2 (x0 + h 2 )2 ] = 0, K4x = hFx (x0 + h, t0 +K3x) = h [ λ1(t0 +K3x) P1 (1− (x0 + h))2 + λ2(t0 +K3x) P2 (x0 + h)2 ] = 0. (35) Then X = x0. (λ1)1 = (λ1)0 + h 6 (K1λ1 + 2K2λ1 + 2K3λ1 +K4λ1) (36) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 11 of 22 Where K1λ1 = hFλ1(x0, t0) = h [ (λ1(t) 2)0 P1 (1− x0)− λ1,0(t)(λ2(t))0 P2 x0(t)− ϕ1 ] = −hϕ1, K2λ1 = hFλ1 ( x0 + h 2 , t0 + K1λ1 2 ) = h [λ1 ( t0 + K1λ1 2 )2 P1 ( 1− ( x0 + h 2 )) − λ1 ( t0 + K1λ1 2 ) λ2 ( t0 + K1λ1 2 ) P2 ( x0 + h 2 ) − ϕ1 ] = −hϕ1, K3λ1 = hFλ1 ( x0 + h 2 , t0 + K2λ1 2 ) = h [λ1 ( t0 + K2λ1 2 )2 P1 ( 1− ( x0 + h 2 )) − λ1 ( t0 + K2λ1 2 ) λ2 ( t0 + K2λ1 2 ) P2 ( x0 + h 2 ) − ϕ1 ] = −hϕ1, K4λ1 = hFλ1 (x0 + h, t0 +K3λ1) = h [ λ1 (t0 +K3λ1) 2 P1 (1− (x0 + h)) − λ1 (t0 +K3λ1)λ2 (t0 +K3λ1) P2 (x0 + h)− ϕ1 ] = −hϕ1. (37) Then λ1,1 = −h2ϕ1 λ2,n+1 = λ2,n + h 6 (K1λ2 + 2K2λ2 + 2K3λ2 +K4λ2) (38) K1λ2 = hFλ2(xn, tn), K2λ2 = hFλ2 ( xn + h 2 , tn + K1λ2 2 ) , K3λ2 = hFλ2 ( xn + h 2 , tn + K2λ2 2 ) , K4λ2 = hFλ2 (xn + h, tn +K3λ2) . (39) Where A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 12 of 22 K1λ2 = hFλ2(x0, t0) = h [ λ1,0(t)(λ2(t))0 P1 (1− x0)− λ2 2,0(t) P2 x0(t) + ϕ1 ] = hϕ2, K2λ2 = hFλ2 ( x0 + h 2 , t0 + K1λ2 2 ) = h [λ1 ( t0 + K1λ2 2 ) λ2 ( t0 + K1λ2 2 ) P1 ( 1− ( x0 + h 2 )) − λ2 ( t0 + K1λ2 2 )2 P2 ( x0 + h 2 ) + ϕ1 ] = hϕ2, K3λ2 = hFλ2 ( x0 + h 2 , t0 + K2λ2 2 ) = h [λ1 ( t0 + K2λ2 2 ) λ2 ( t0 + K2λ2 2 ) P1 ( 1− ( x0 + h 2 )) − λ2 ( t0 + K2λ2 2 )2 P2 ( x0 + h 2 ) + ϕ1 ] = hϕ2, K4λ2 = hFλ2 (x0 + h, t0 +K3λ2) = h [ λ1 (t0 +K3λ2)λ2 (t0 +K3λ2) P1 (1− (x0 + h)) − λ2 (t0 +K3λ2) 2 P2 (x0 + h) + ϕ1 ] = hϕ2. (40) Then λ2,1 = h2ϕ2 u1,1 = λ1,1 P1 (1− x1) = −h2ϕ1 P1 (1− x0), u2,1 = −λ2 P2 x = −h2ϕ2 P2 x, c1,1 = p1 2 u21,1 = p1 2 ( h2ϕ1 P1 (1− x0)) 2, c2,1 = p2 2 u22,1 = p2 2 ( h2ϕ2 P2 x)2. (41) At n=1 x2 = x1 + h 6 (K1x + 2K2x + 2K3x +K4x) (42) Where A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 13 of 22 K1x = hFx(x1, t1) = h [ (λ1(t))1 P1 (1− x1(t)) 2 + (λ2(t))1 P2 x21 ] , K2x = hFx ( x1 + h 2 , t1 + K1x 2 ) = h [ λ1 ( t1 + K1x 2 ) P1 ( 1− ( x1 + h 2 ))2 + λ2 ( t1 + K1x 2 ) P2 ( x1 + h 2 )2 ] , K3x = hFx ( x1 + h 2 , t1 + K2x 2 ) = h [ λ1 ( t1 + K2x 2 ) P1 ( 1− ( x1 + h 2 ))2 + λ2 ( t1 + K2x 2 ) P2 ( x1 + h 2 )2 ] , K4x = hFx (x1 + h, t1 +K3x) = h [ λ1 (t1 +K3x) P1 (1− (x1 + h))2 + λ2 (t1 +K3x) P2 (x1 + h)2 ] . (43) Then x2 = x0+ h4 6 [ ϕ1 p1 ( (1−x0) 2+4(1−x0− h 2 )2+(1−x0−h)2 ) −ϕ2 p2 ( (x0) 2+4(x0+ h 2 )2+(x0+h)2 )] . (λ1)2 = λ1,1 + h 6 (K1λ1 + 2K2λ1 + 2K3λ1 +K4λ1) (44) Where K1λ1 = hFλ1(x1, t1) = h [ (λ2 1(t))1 P1 (1− x1)− (λ1(t))1(λ2(t))1 P2 x1(t)− ϕ1 ] , K2λ1 = hFλ1 ( x1 + h 2 , t1 + K1λ1 2 ) = h [ λ1(t1 + K1λ1 2 )2 P1 (1− (x1 + h 2 )) − λ1(t1 + K1λ1 2 )λ2(t1 + K1λ1 2 ) P2 (x1 + h 2 )− ϕ1 ] , K3λ1 = hFλ1 ( x1 + h 2 , t1 + K2λ1 2 ) = h [ λ1(t1 + K2λ1 2 )2 P1 (1− (x1 + h 2 )) − λ1(t1 + K2λ1 2 )λ2(t1 + K2λ1 2 ) P2 (x1 + h 2 )− ϕ1 ] , K4λ1 = hFλ1 (x1 + h, t1 +K3x) = h [ λ1(t1 +K3λ1) 2 P1 (1− (x1 + h)) − λ1(t1 +K3λ1)λ2(t1 +K3λ1) P2 (x1 + h)− ϕ1 ] . (45) λ1,2 = −h2ϕ1 + h6 6 [ (ϕ1) 2 P1 ( (1− x0) + 4(1− x0 − h 2 ) + (1− x0 − h) ) − ϕ1ϕ2 P2 ( x0 + 4(x0 + h 2 ) + (x0 + h) )] + ϕ1h 2 (46) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 14 of 22 K1λ1 = hFλ1(x1, t1) = h [ λ1,1(t) P1 (1− x1)− λ1,1(t)λ2,1(t) P2 x1(t)− ϕ1 ] , K2λ1 = hFλ1 ( x1 + h 2 , t1 + K1λ1 2 ) = h [ λ1 ( t1 + K1λ1 2 )2 P1 ( 1− (x1 + h 2 ) ) − λ1 ( t1 + K1λ1 2 ) λ2 ( t1 + K1λ1 2 ) P2 (x1 + h 2 )− ϕ1 ] , K3λ1 = hFλ1 ( x1 + h 2 , t1 + K2λ1 2 ) = h [ λ1 ( t1 + K2λ1 2 )2 P1 ( 1− (x1 + h 2 ) ) − λ1 ( t1 + K2λ1 2 ) λ2 ( t1 + K2λ1 2 ) P2 (x1 + h 2 )− ϕ1 ] , K4λ1 = hFλ1 (x1 + h, t1 +K3λ1) = h [ λ1 ( t1 +K3λ1 )2 P1 (1− (x1 + h)) − λ1 ( t1 +K3λ1 ) λ2 ( t1 +K3λ1 ) P2 (x1 + h)− ϕ1 ] . (47) Similarly, λ2,2 = h2ϕ2 − h6 6 [ ϕ1ϕ2 P1 ( (1− x0) + 4(1− x0 − h 2 ) + (1− x0 − h) ) + (ϕ2) 2 P2 ( x0 + 4(x0 + h 2 ) + (x0 + h) )] + ϕ2h 2. (48) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 15 of 22 u1,2 = λ1,2 P1 (1− x2) = 1 P1 [( 1− ( x0 + h4 6 [ ϕ1 p1 ( (1− x0) 2 + 4(1− x0 − h 2 )2 + (1− x0 − h)2 ) − ϕ2 p2 ( x20 + 4(x0 + h 2 )2 + (x0 + h)2 )]))] · [ h2ϕ1 + h6 6 ( (ϕ1) 2 P1 ( (1− x0) + 4(1− x0 − h 2 ) + (1− x0 − h) ) − ϕ1ϕ2 P2 ( x0 + 4(x0 + h 2 ) + (x0 + h) )) + ϕ1h 2 ] , u2,2 = −λ2,2 P2 x2 = − 1 P2 [ − h2ϕ2 − h6 6 ( ϕ1ϕ2 P1 ( (1− x0) + 4(1− x0 − h 2 ) + (1− x0 − h) ) + (ϕ2) 2 P2 ( x0 + 4(x0 + h 2 ) + (x0 + h) )) + ϕ2h 2 ] · [ x1 + h 6 (K1x + 2K2x + 2K3x +K4x) ] , (49) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 16 of 22 c1,2 = p1 2 (u21)2 = p1 2 [ 1 P1 [ 1− ( x0 + h4 6 [ ϕ1 p1 ( (1− x0) 2 + 4(1− x0 − h 2 )2 + (1− x0 − h)2 ) − ϕ2 p2 ( x20 + 4(x0 + h 2 )2 + (x0 + h)2 )])] · [ h2ϕ1 + h6 6 ( (ϕ1) 2 P1 ( (1− x0) + 4(1− x0 − h 2 ) + (1− x0 − h) ) − ϕ1ϕ2 P2 ( x0 + 4(x0 + h 2 ) + (x0 + h) )) + ϕ1h 2 ]]2 , c2,2 = p2 2 (u22)2 = p2 2 [ − 1 P2 [ − h2ϕ2 − h6 6 ( ϕ1ϕ2 P1 ( (1− x0) + 4(1− x0 − h 2 ) + (1− x0 − h) ) + (ϕ2) 2 P2 ( x0 + 4(x0 + h 2 ) + (x0 + h) )) + ϕ2h 2 ] · [ x1 + h 6 (K1x + 2K2x + 2K3x +K4x) ]]2 (50) Suppose that Firm 2 is already present in the market, with its market share represented by 1 − x(t). Firm 1 enters the market to compete with Firm 2, and its market share at time t is represented by x(t). We assume that at the beginning of the competition, the initial market share of Firm 1 is x0 = 0. The parameters are set as follows: ϕ1 = 0.1, ϕ2 = 0.3, P1 = 0.25, P2 = 0.5, the time interval t ∈ [0, 1.5], and the step size h = 0.5. Then we have a comparison between the two firms. At n=0, x0 = 0, 1− x0 = 1, u1,0 = 0, u2,0 = 0, c1,0 = 0, c2,0 = 0, λ1,0 = 0, λ2,0 = 0. (51) This indicates that, initially, Firm 2 was the only firm in the market and did not undertake any advertising efforts. At n=1 x1 = 0, 1− x1 = 1, u1,1 = 0.1, u2,1 = 0, c1,1 = 0.00125, c2,1 = 0, λ1,1 = 0.25, λ2,1 = −0.075. (52) A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 17 of 22 This indicates that initially, Firm 1 entered the market without holding any market share. Consequently, it launched advertising campaigns to increase its market share, while Firm 2 did not engage in any advertising. As a result, the number of customers for Firm 1 increased. At n=2, x2 = 0.00989583, 1− x2 = 0.99010417, u1,2 = 0.18533512, u2,2 = −0.0001182685, c1,2 = 0.004293638, c2,2 = 3.496865× 10−9, λ1,2 = 0.04679687, λ2,2 = −0.005976. (53) This indicates that initially, Firm 1’s market share increased, while Firm 2 began to lose market share and experienced a decrease in sales. As a result of these losses, Firm 2 obtained a loan to finance its advertising campaigns. At n=3, x3 = 0.0362395, 1− x3 = 0.96376, u1,3 = 0.283165810, u2,3 = 0.010028594, c1,3 = 0.005976423, c2,3 = 0.00446470, λ1,3 = 0.0734531, λ2,3 = −0.082457305. (54) Now, we present a comparison of the obtained numerical solutions for the two firms in the following figures: Figure 1 shows the market shares x and 1 − x, representing Firm 1 and Firm 2, respectively, at time t. It is observed that Firm 1’s market share increases over time, while Firm 2’s market share decreases. Figure 2 illustrates the controls u1 and u2, representing the advertising efforts of Firm 1 and Firm 2, respectively. Initially, Firm 1 increased its advertising efforts to attract more customers, while Firm 2 had no advertising. Over time, Firm 2 gradually increased its advertising to reduce Firm 1’s market share. Figure 3 presents the advertising cost functions for the two firms, c1 for Firm 1 and c2 for Firm 2. Firm 1’s cost is noticeably higher due to its stronger advertising intensity, whereas Firm 2 incurs minimal costs during the early stages of the competition. Figure 4 depicts the costate variables λ1 and λ2 for Firm 1 and Firm 2, respectively. For Firm 1, the costate variable starts at zero since it is new to the market and initially has no impact. It then rises as Firm 1 gains influence and faces competition, and later decreases as it approaches the optimal strategy, indicating a balanced state. Conversely, λ2 initially starts at zero and then decreases to negative values as Firm 2 loses market share due to Firm 1’s entry. The negative values indicate that maintaining its market share becomes increasingly costly, highlighting the competitive pressure imposed by Firm 1. A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 18 of 22 1 2 3 4 t 0.2 0.4 0.6 0.8 1.0 x Figure 1: Market shares of Firm 1 and Firm 2 at time t. 0.5 1.0 1.5 2.0 2.5 3.0 t 0.05 0.10 0.15 0.20 0.25 U Figure 2: The advertising efforts of Firm 1 and Firm 2 at time t. 0.5 1.0 1.5 2.0 2.5 3.0 t 0.002 0.004 0.006 0.008 0.010 C Figure 3: The advertising cost functions of Firm 1 and Firm 2 at the time t 0.5 1.0 1.5 2.0 2.5 3.0 -0.05 0.05 0.10 0.15 0.20 0.25 Figure 4: The costate variables of Firm 1 and Firm 2 at time t. A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 19 of 22 4. Stability Analysis of the Fourth-Order Runge-Kutta Method Using the Jacobian Matrix In a system of differential equations, stability analysis is performed to determine whether small perturbations around equilibrium points cause the system to return to equilibrium (stable) or to diverge away (unstable). For a nonlinear system: dx dt = λ1 P1 (1− x)2 + λ2 P2 x2, dλ1(t) dt = λ2 1 P1 (1− x)− λ1λ2(t) P2 x− ϕ1, dλ2(t) dt = λ1λ2(t) P1 (1− x)− λ2 2 P2 x+ ϕ2. (55) Equilibrium points (or steady states) of the system are obtained by solving: f1 = λ1 P1 (1− x)2 + λ2 P2 x2 = 0, f2 = λ2 1 P1 (1− x)− λ1λ2(t) P2 x− ϕ1 = 0, f3 = λ1λ2(t) P1 (1− x)− λ2 2 P2 x+ ϕ2 = 0. (56) After solving the system (56) simultaneously, we obtain the equilibrium point (x∗, λ∗ 1, λ ∗ 2). Next, we construct the Jacobian matrix, which is a square matrix of the first-order partial derivatives of the system’s functions with respect to the state variables. The Jacobian matrix J is defined as: J =  ∂f1 ∂x ∂f1 ∂λ1 ∂f1 ∂λ2 ∂f2 ∂x ∂f2 ∂λ1 ∂f2 ∂λ2 ∂f3 ∂x ∂f3 ∂λ1 ∂f3 ∂λ2  , Then the jacobian matrix evaluation, J =  −2λ∗ 1 p1 (1− x∗) + 2λ∗ 2 p2 x (1−x∗)2 p1 (x∗)2 p2 −(λ∗ 1) 2 p1 − λ1λ∗ 2 p2 2λ∗ 1(1−x∗) p1 − λ2(t) P2 x −λ∗ 1x ∗ p2 −λ∗ 1λ ∗ 2 p1 − (λ∗ 2) 2 p2 λ∗ 2(1−x∗) p1 λ∗ 1(1−x∗) p1 − 2λ∗ 2x ∗ p2  . To analyze the stability at an equilibrium point (x∗, λ∗ 1, λ ∗ 2), the Jacobian matrix J(x∗, λ∗ 1, λ ∗ 2) is evaluated at the equilibrium. The stability of the system depends on the eigenvalues µi of J(x ∗, λ∗ 1, λ ∗ 2). Stable Equilibrium: All eigenvalues satisfy Re(µi) < 0, implying local asymptotic stability. Unstable Equilibrium: At least one eigenvalue satisfies Re(µi) > 0, implying instabil- ity. A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 20 of 22 Marginally Stable Equilibrium: All eigenvalues satisfy Re(µi) = 0; stability cannot be fully determined and the system may be marginally stable or exhibit oscillatory behavior. To study the stability of the system (55), we first determine the equilibrium points by solving the system (56). We obtain the following equilibrium points: (x∗, λ∗ 1, λ ∗ 2)1 = (1.85786, 0− 0.617178i, 0 + 0.263177i), (x∗, λ∗ 1, λ ∗ 2)2 = (1.85786, 0 + 0.617178i, 0− 0.263177i), (x∗, λ∗ 1, λ ∗ 2)3 = (−0.916188, 0 + 0.0640229i, 0− 0.560109i), (x∗, λ∗ 1, λ ∗ 2)4 = (−0.916188, 0− 0.0640229i, 0 + 0.560109i), (x∗, λ∗ 1, λ ∗ 2)5 = (0.391662, 0.100038,−0.482684), (x∗, λ∗ 1, λ ∗ 2)6 = (0.391662,−0.100038, 0.482684), Then the Jacobian matrix evaluation, J =  −2λ∗ 1 0.25 (1− x∗) + 2λ∗ 2 0.5 x − (1−x∗)2 0.25 (x∗)2 0.5 −(λ∗ 1) 2 0.25 − λ∗ 1λ ∗ 2 0.5 2λ∗ 1(1−x∗) 0.25 −λ∗ 1x ∗ 0.5 −λ∗ 1λ ∗ 2 0.25 − (λ∗ 2) 2 0.5 λ∗ 2(1−x∗) 0.25 λ∗ 1(1−x∗) 0.25 − 2λ∗ 2x ∗ 0.5  . The eigenvalues µi of the Jacobian matrix are computed, J(x∗, λ∗ 1, λ ∗ 2)1 =(−2.23227− 4.23449i,−2.23227− 0.97901i,−6.85728× 10−15 − 2.279837i), J(x∗, λ∗ 1, λ ∗ 2)2 =(−2.23227 + 4.23449i,−2.23227 + 0.97901i,−6.85728× 10−15 + 2.279837i), J(x∗, λ∗ 1, λ ∗ 2)3 =(−1.04710− 2.79990i,−1.04710− 0.79213i, 8.88260× 10−16 − 1.07122i), J(x∗, λ∗ 1, λ ∗ 2)4 =(−1.04710 + 2.79990i,−1.04710 + 0.79213i, 8.88260× 10−16 + 1.07122i), J(x∗, λ∗ 1, λ ∗ 2)5 =(−1.24304,−0.62152,−0.51276), J(x∗, λ∗ 1, λ ∗ 2)6 =(−1.24304,−1.24304,−0.62152). From the computed eigenvalues of the Jacobian matrix, all eigenvalues have negative real parts, indicating that the equilibrium point of the system is locally asymptotically stable. This means that small perturbations in the market shares or advertising efforts of either firm will decay over time, and the system will return to the equilibrium state. Economically, this implies that both Firm 1 and Firm 2 can maintain stable market shares under the given dynamics, and neither firm will experience unbounded fluctuations in influence or costs. 5. Conclusion In this paper, we discussed the numerical solution of an open-loop Nash differential game using the fourth-order Runge-Kutta method. This approach allowed us to analyze competition between firms in a dynamic market environment. We studied the convergence A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 21 of 22 of the Runge-Kutta method to ensure the reliability of the numerical solution, and we examined the stability of the solution to confirm that numerical errors do not amplify over time. We also demonstrated that increasing advertising campaigns leads to higher market participation for the firms. Finally, illustrative figures were presented to demonstrate the effect of advertising campaigns on each firm’s market share. In future work, we plan to expand this study in several ways. First, we aim to incorporate stochastic elements into the model to analyze how uncertainty influences market behavior and firm strategies. Second, we plan to include more firms, allowing multiple companies to compete simultaneously, which increases the complexity of interactions and may reveal new patterns. Third, we aim to allow strategies to adapt based on the current state of the system, using feedback Nash equilibria, making the model more flexible and closer to real-world competition. References [1] John Von Neumann and Oskar Morgenstern. Theory of Games and Economic Be- havior. Princeton University Press, 1944. [2] Fred Brauer and John A. Nohel. Differential Equations and Their Applications. Springer, 1989. [3] Mojtaba Dehghan Banadaki and Hamidreza Navidi. Numerical solution of open-loop nash differential games based on the legendre tau method. Games, 11:28, 2020. [4] Dmitry R. Kuvshinov and Sergei I. Osipov. Numerical construction of stackelberg solutions in a linear positional differential game based on the method of polyhedra. Automation and Remote Control, 79:479–491, 2018. [5] A.A. Megahed and S.E. Mohamed. On the solution of differential games by picard method. Applied Mathematics and Computation, 184:432–440, 2007. [6] A.A. Megahed and A.M. Kassem. A zero-sum differential game model for the spread of coronavirus. Chaos, Solitons & Fractals, 140:110–117, 2020. [7] E. Youness and A.M.A. El-Sayed. A stackelberg differential game approach to combat terrorism. Nonlinear Dynamics, 63:671–682, 2011. [8] A.A. Hemeda. On iterative methods for nonlinear integro-differential equations. Ap- plied Mathematics and Computation, 219:7625–7635, 2013. [9] A. A. Megahed, A. A. Hemeda, and H. F. A. Madkour. Approximate solutions for nash differential games. Journal of Mathematics, 2023:7374882, 2023. [10] A.M. Kassem and A.A. Megahed. A differential game approach to counter-terrorism. Applied Mathematics and Computation, 190:607–618, 2007. [11] E. Youness. Nash-collaborative solutions for differential games. Applied Mathematics and Computation, 138:289–300, 2003. [12] E. Youness and A.A. Megahed. On nash equilibrium strategies of fuzzy differential games. Fuzzy Sets and Systems, 113:439–445, 2000. [13] H. Sun and J. Yong. Open-loop and closed-loop nash equilibria in linear-quadratic stochastic differential games. SIAM Journal on Control and Optimization, 41:737– 756, 2002. A. A. Megahed et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6007 22 of 22 [14] J.C. Engwerda. Open-loop nash equilibria in lq-games. Journal of Economic Dynam- ics and Control, 22:29–50, 1998. [15] E. Youness. A fuzzy min-max differential game problem. Fuzzy Sets and Systems, 122:465–474, 2001. [16] E. Youness. Large scale fuzzy differential games. Applied Mathematics and Compu- tation, 142:225–234, 2003. [17] A.A. Megahed. A min-max differential game applied to terrorism. Applied Mathe- matics and Computation, 198:446–455, 2008. [18] Ebrahim A. Youness et al. Min-max differential game with partial differential equa- tion. AIMS Mathematics, 7:13777–13789, 2022. [19] Kendall E. Atkinson. Numerical Methods for Ordinary Differential Equations. Wiley, 1989. [20] Ebrahim Youness and Abd El-Moneim Megahed. A study on fuzzy differential game. Le Matematiche, 56:97–107, 2001.