Electronic Journal of Differential Equations, Vol. 2023 (2023), No. 59, pp. 1–22. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu DOI: 10.58997/ejde.2023.59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL WITH SOCIAL MOBILITY, RENEWABLES AND NONRENEWABLES MEIR SHILLOR, THANAA ALI KADHIM Abstract. We expand the HANDY (Human And Nature DYnamics) model for the socioeconomic dynamics of a large stratified society. The basic model was introduced in Motesharrei (Dissertation 2014) and Motesharrei et al. (2016). It is a nonlinear system of ODEs for a ‘simple society’ of Elites, Workers, Wealth, and Natural Resources. Following Ali Kadhim (Dissertation 2021), we add social mobility between the classes and split natural resources into renewables and nonrenewables. We establish the existence, bounded- ness and positivity of the solutions, and investigate the stability of the steady states. The model admits stable steady states, and there is numerical evidence of stable periodic solutions and limit cycles. Simulations depict the different qualitative types of model behavior: convergence to steady states, periodic oscillations or collapse. 1. Introduction The HANDY – Human And Nature DYnamics – was a ground breaking sus- tainability focused mathematical model for complex societies, constructed in Mote- sharrei’s dissertation [5]. The original model consisted of four ordinary differential equations (ODEs), representing the growth rates of two populations, their use of natural resources and of their wealth. The model was simulated to study its be- havior in [5] and in Motesharrei et al. [6, 7]. Computer simulations of the HANDY model showed that economic stratification and rapid depletion of natural resources were among the main reasons for societal collapse. The simulation results indicated that an irreversible collapse could be avoided if the society adopted sustainable be- haviors to reduce the depletion of natural resources and moved towards an equitable distribution of resources. In these publications, the authors provided many histor- ical examples and studies of collapse of earlier civilizations, looking for possible underlying processes that were common to all, instead of trying to study each so- ciety separately. Indeed, the main interest in this ‘simple model’ and its computer simulations was to find common causes for the collapse of the civilizations, in addi- tion to issues specific to each one. It was found that the common causes in all the cases were both over-depletion of resources and inequality. This finding of common 2020 Mathematics Subject Classification. 91D10, 37N99, 34C60. Key words and phrases. HANDY model; social mobility’ renewables and nonrenewables; socioeconomic dynamical system; ODEs analysis. ©2023. This work is licensed under a CC BY 4.0 license. Submitted January 17, 2023. Published September 15, 2023. 1 2 M. SHILLOR, T. A. KADHIM EJDE-2023/59 causes for the collapse of the different societies is what makes mathematical models so important and effective in clarifying some general types of behavior that are hard to see when studying each society separately. Related previous work is as follows. In Brendan and Taylor [2], the authors introduced a simpler model to study the collapse of the society in Easter Island. They used the Ricardo-Malthus model with two ODE’s, one for renewable resources that included a logistic term with saturation and a harvesting term. The popu- lation growth rate was described by three terms: the birth rate, the death rate, and a fertility factor using the Malthusian term. The latter claims that higher income rate stimulates population growth, leading to the depletion of resources and eventually to famine. On the other hand, it was found that too low wages may cause a population collapse caused by conflict. Various civilizations have gone through such collapse, including the western Roman empire and the Mayan and the Mesopotamian civilizations. A similar model for the societal dynamics was considered by Grammaticos et al. in [3], consisting of three ODEs representing the population, resources and reserves. They used a special factor for the amount of reserve per capita, and proposed that the population was flourishing when the accumulated reserve was abundant (above the threshold) and died out when the reserves were insufficient. Furthermore, they found in their simulations patterns, similar to those in the HANDY model, of nature’s slow growth and rapid decline, known as the Seneca effect. This work continues the study of such models, and is based on the dissertation [1], where a few different versions of the model were constructed, analyzed and computer simulated. Here, we have a three-fold aim: First, we expand the HANDY model by allowing for social mobility between the workers and the elite, and by splitting the resources into renewable energy/wealth sources (wood, solar, wind), and nonrenewable sources (oil, gas, coal). Secondly, we perform mathematical analysis of the model, proving the existence, boundedness and positivity of the solutions, as well as a study of the system steady states, following [1]. Thirdly, we present simulations of typical types of model predictions, and initial study of the dependence of the solutions on the consumption rates of the resources and on inequity in consumption. The computer simulations indicate that the system has four qualitatively dif- ferent types of behavior: monotone approach to the steady state, Figure 1 (a) (where the approach is monotone after one initial oscillation); damped oscillations approach, Figure 1 (b), periodic oscillations, Figure 1 (c), and system collapse, Figure 1 (d). The different cases depend on the different sets of initial conditions and the values of the various system parameters, and provide insight into possible behaviors of such models and hence, to the potential behavior of complex societies. It seems that such extensions of the HANDY model may be useful in investigating various scenarios and the associated conditions for flourishing or collapse of com- plex societies. It is very likely that the system is rich enough to exhibit chaotic solutions, however, our simulations did not find any. Moreover, we do not have a mathematical proof of the existence of periodic solutions, yet. Section 3 establishes the existence, positivity and boundedness of the solutions of the model. Similar properties of the HANDY model have been established in [1], and the analysis here is tailored to the HANDY-SMRN model. Section 4, studies the system’s steady states and their stability. It is found that the steady EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 3 states depend essentially on the ratio η = w/wth of the wealth w to the wealth’s threshold wth. Because of the structure of the system coefficients, there are three regimes, which lead to three families of steady states: Case (i), 0 ≤ η ≤ 1/κ1, which includes the two steady states with collapse, namely (0, 0, 0, 0, 0, ) and (0, 0, λ, 0, 0). The first of which is found to be unstable while the second is stable and attracting. Case (ii), 1/κ1 ≤ η < 1, where the elite have enough wealth; and Case (iii), 1 < η, where there is enough wealth (which includes food reserves) for both populations, to avoid collapse. We note that the system has 15 parameters, and five initial conditions, so a complete numerical study is out of reach. Therefore, we present only the four qualitatively different types of solutions; simulations with different δn and δr; and with different κ1. We note that a more thorough computer study of the effects of different rates of depletion of renewables and nonrenewables may be of interest. The computer simulations with different values of the depletion of nonrenew- ables, with two different values of δn show that a more moderate use of nonrenew- ables leads to longer periods of reaching the steady states but also to substantially higher steady state populations. Concerning the renewable natural resources, with depletion factor δr, we found that there is a range of values for which the solutions are either oscillatory, or oscillate for very long periods of time and then quickly converge to the steady states. Finally, simulations with different consumption rate constant κ1 show that for κ1 = 1, 5, 10 the approach to the steady states is increasingly more complex with more large oscillations, in which the populations almost die out. Note: It is an important observation that even in the cases where, eventually, the solutions converge to similar steady states, the ways these are reached may be very different, with more or less societal upheaval on the way. This indicates that just studying the steady states of the model and their stability may provide only a very partial and limited description of the society’s dynamics. Clearly, there is considerable interest to find out if the model can fit some of the historical collapses of societies which essentially consisted of commoners and nobility. The rest of the work is organized as follows. The model is presented in Section 2, where the underlying assumptions are discussed with some detail, taken from the original HANDY model and [1]. Then, we discuss our expansion of the model by allowing for mobility between the two classes, and splitting the natural resources into renewables and nonrenewables. We also introduce the ratio of wealth to the threshold wealth, η, which plays an important role in the analysis. Section 3 deals with the analysis of the model. It establishes the boundedness and positivity of the solutions, Proposition 3.2. These lead to the existence of the solutions on every finite time interval, Theorem 3.3. Then, Section 4 studies the system steady states and establishes the conditions for their stability. Section 5 presents the results of model simulations. It depict the four main types of behavior that the model predicts. Also it presents a short study of the effects of increasing or decreasing the renewables and nonrenewables rate constants, and the dependence on the consumption rate constant. Section 6 poses some possible unresolved questions of interest. 4 M. SHILLOR, T. A. KADHIM EJDE-2023/59 2. HANDY-SMRN model The extended model studied here, is a version of the model constructed in Ali Kadhim [1] (Dissertation 2021; Chs. 5 and 6) and consists of five coupled nonlinear ODEs: two equations model the populations growth rates of the Commoners or workers, xc(t), and the Elites or the rich, xe(t); one describes the rate of growth or depletion of renewable natural resources, yr(t) (such as wood, hydro-power, or solar energy), and one for the depletion of non-renewables (such as coal, oil or gas) yn(t), and the rate of growth of wealth w(t), which includes food surpluses. The time t is measured in ‘time units,’ (in real societies this would be in years). The novelties in this model, as noted in the introduction, are the addition of social mobility between the rich and the workers, and the splitting of the natural resources into renewables and nonrenewables. The complete HANDY-SMRN model, the HANDY model with Social Mobility, Renewables and Nonrenewables, is the following: Problem 2.1. Find five functions (xc(t), xe(t), yr(t), yn(t), w(t)), defined on [0, T ], 0 < T <∞, such that dxc dt = (βc − αc)xc + γexe − γcxc, (2.1) dxe dt = (βe − αe)xe − γexe + γcxc, (2.2) dyr dt = γyr(λ− yr)− δrxcyr, (2.3) dyn dt = −δnxcyn, (2.4) dw dt = δwxc(yr + yn)− Cc − Ce. (2.5) Here, the consumption rate functions Cc and Ce are given in (2.9), the wealth threshold wth is given in (2.7), and the death rates αc and αe are given in (2.10). The initial conditions are, xc(0) = xc0, xe(0) = xe0, yr(0) = yr0, yn(0) = yn0 w(0) = w0. (2.6) We assume that xc0 > 0, xe0 ≥ 0, yr0 > 0, yn0 ≥ 0, w0 ≥ 0, are prescribed. The rest of the notation is as follows. The respective birth rates are βc, βe, assumed to be constant; γ is the nature’s regeneration factor, λ its saturation level or the carrying capacity of the renewable resources; δr and δn are the respective renewable and non-renewable resources depletion rate constants; δw is the conversion factor between the resources depletion and wealth generation. Next, γe is the rate constant in which the rich become workers or commoners (going bankrupt), and γc is the rate constant in which the workers ‘make it’ and become wealthy. All these may be constants or given non-negative functions. We note that (2.3) is of the logistic type, and allows for regeneration of the renewable resources up to the carrying capacity λ, while equation (2.4) for the nonrenewable does not allow for regeneration. Before we describe Cc, Ce, the total consumption rates, we introduce the wealth threshold wth = ρ(xc + κ0xe + ε), (2.7) EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 5 where ε > 0 is a very small number that prevents the threshold wth from becoming zero, since it appears in the denominators in (2.9) and (2.10). We next introduce the variable which plays an important role in the analysis of the steady states, η = η(t) = w(t) wth(t) . (2.8) The total consumption rates are given by Cc = min(1, η)sxc, Ce = min(1, η)κ1sxe. (2.9) Finally, the death rates αc, αe, are given by, αc = αm + max(0, 1−min(1, η))(αM − αm), αe = αm + max(0, 1− κ1 min(1, η))(αM − αm). (2.10) We refer to [5] and to [1] for a more detailed description of the model parameters and coefficients. We note that in the original HANDY model [5], the choice was y = yr + yn, κ0 = κ1 = κ, γc = γe = 0, δr = δw = δ, δn = 0. 3. Model analysis This section provides the analysis of the model, starting with establishing the boundedness, positivity and then the existence of the solutions to Problem 2.1. Then, we study the steady states of the system and establish their stability. The existence of the solution is based on the following theorem that can be found in Kuttler [4], and for the benefit of the reader, its statement is included. Theorem 3.1. Let F : [0, T ] × Rn → Rn be continuous and suppose that there exists L > 0 such that for all λ̂ ∈ (0, 1), if x′ = λ̂F (x, t) , x(0) = x0, (3.1) for all t ∈ [0, T ], implies that ‖x‖ < L. Then, there exists a solution to the system of ODEs, x′ = F(x, t), x(0) = x0. (3.2) for t ∈ [0, T ]. We note that the theorem does not guarantee the uniqueness of the solution. To show uniqueness, a more refined argument is needed, which we leave open. 3.1. Positivity and a priori estimates. Let z = (xc, xe, yr, yn, w) and F : R5 → R5 be given by F(z) =  βcxc − αcxc + γexe − γcxc βexe − αexe − γexe + γcxc γyr(λ− yr)− δrxcyr −δnxcyn δwxc(yr + yn)− Cc − Ce,  =  F1 F2 F3 F4 F5  . (3.3) Here, Cc and Ce are given in (2.9), wth is given in (2.7), and αc and αe are given in (2.10). We assume, following [5, 1], that 0 < αm < βe ≤ βc ≤ αM < 1, and then it follows that αM − αm βe − αm ≥ αM − αm βc − αm ≥ 1. (3.4) Everywhere below, we assume that κ0, κ1 ≥ 1. 6 M. SHILLOR, T. A. KADHIM EJDE-2023/59 First, we write the system as z′ = F (z) , z(0) = z0 = ( xc0, xe0, yr0, yn0, w0 ) . (3.5) Here and below, the prime indicates the time derivative, and each component of z0 is assumed to be positive. Below, we show that the results hold true for xe0 = 0, yn = 0 and w0 = 0, as well. We next establish the necessary a-priori estimates on the system solutions. Let 0 < T < ∞, be fixed, and then [0, T ] is an arbitrary finite interval. We begin by showing that starting with positive initial conditions, all the possible solutions are nonnegative and bounded. This leads to estimates that are valid for [0, T ] and allow us to prove that a solution exists on this interval. Proposition 3.2. Assume that xc0 > 0, yr0 > 0, yn0 > 0 and xe0, w0 are non- negative. If z = (xc, xe, yr, yn, w) is any solution to the system (3.5), then all the components of z are nonnegative. Moreover, for 0 ≤ t ≤ T , the following estimates hold: 0 ≤ xc(t) ≤ (xc0 + xe0)eβcT , 0 ≤ xe(t) ≤ (xc0 + xe0)eβcT , 0 < yr(t) ≤ yr0eγλT , 0 < yn ≤ yn0, 0 ≤ w(t) ≤ (w0 +M) +Me(γλ+βc)T . (3.6) Here, M is a positive constant (see below) that depends only on the problem parameters. Proof. It follows from the first two equations, by multiplying the first with xe and the second with xc, adding the expressions and rearranging, that (xexc) ′ = ((βc − αc) + (βe − αe))xcxe − (γc + γe)xexc + γcx 2 c + γex 2 e. Since initially xe(0)xc(0) > 0, by continuity, there is t0 > 0 such that xe(t)xc(t) > 0 for 0 ≤ t < t0. Then, if xe(t0)xc(t0) = 0, we obtain( xexc )′ (t0) = γcx 2 c + γex 2 e ≥ 0. If either xe(t0) > 0 or xc(t0) > 0, then (xexc) ′(t0) > 0, so, by continuity, the func- tion (xexc)(t) is increasing at t0, which contradicts the assumption that xe(t0)xc(t0) = 0 for the first time. Now, suppose that both vanish at t0 > 0. Then, for t ≥ t0, we consider the system x′c(t) = (βc − αc)xc + γexe − γcxc, xc(t0) = 0, x′e(t) = (βe − αe)xe − γexe + γcxc, xe(t0) = 0, which has the unique trivial solution solution xc(t) = xe(t) = 0 for t ≥ t0, so extinction has occurred at t0. It follows that xexc ≥ 0 on [0, T ], for each 0 < T and so both functions are nonnegative, as long as they exist. Next, to obtain upper bounds on xc and xe, we add the first two equations and obtain (xc + xe) ′ = (βc − αc)xc + (βe − αe)xe ≤ βc(xc + xe). Therefore, xc(t) + xe(t) ≤ (xc0 + xe0)eβct. EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 7 Since xc and xe are nonnegative, each one is bounded and this establishes the first two estimates in (3.6). We show now that yr is actually positive. We note that y′r = γyr(λ− yr)− δrxcyr. Since xc ≥ 0 and γ, δr > 0, and since y′r ≤ γλyr, we obtain yr(t) ≤ yr0eγλt. Then, since y′r ≥ −(γyr + δrxc)yr, it follows that y′r ≥ −(γyr0e γλt + δrxcoe βct)yr. Thus, yr(t) ≥ yr0 exp (−(γyr0e γλt + δr(xco + xe0)eβct)) > 0. This proves both sides of the estimate for yr in (3.6). Finally, we address the last inequality for w. We have w′ = δwxcyr − (min (1, η) sxc + min (1, η)κ1sxe) Since w0 > 0, assume that there is 0 < t0 such that w(t) > 0 for 0 ≤ t < t0, and if w(t0) = 0 then η(t0) = 0, and so w′(t0) = δwxcyr ≥ 0. But y > 0 so if xc(t0) = 0, then, as discussed above, both xc and xe vanish identically for t0 ≤ t, which shows that w(t) = 0 for t0 ≤ t, as well. Furthermore, it follows from the equation that w′ ≤ δwxcyr, and using the estimates for yr and xc, we obtain w′ ≤ δwyr0(xc0 + xe0)e(γλ+βc)t. Integration over [0, t] yields w(t) ≤ (w0 +M) +Me(γλ+βc)t, where M = δwyr0(xc0 + xe0) (γλ+ βc) . We conclude that w is nonnegative and bounded. This completes the proof. � These estimates allow us to establish the existence of solutions of the system by using Theorem 3.1. Indeed, the a-priori estimates, established in Proposition 3.2, for possible solutions of Model 2.1, allow us to apply Theorem 3.1, hence we obtain the existence of the solutions to the HANDY-SMRN Model. Moreover, a careful checking of the proof of Proposition 3.2 shows that the results hold true even when x0e = 0 and w0 = 0, by allowing x0e → 0 and w0 → 0. We summarize these results in our main existence theorem. Theorem 3.3. Assume that xc0 and yr0 are positive and xe0, yn0 and w0 are non- negative. Then, Model 2.1 has solutions on every finite time interval [0, T ]. More- over, the solutions satisfy estimates (3.6). As was mentioned above, Theorem 3.1 does not require the right-hand side F to be Lipschitz continuous, but then, Theorem 3.3 does not guarantee the uniqueness of the solution. That has to be shown separately, and is left open in this work. 8 M. SHILLOR, T. A. KADHIM EJDE-2023/59 4. Steady states and stability This section studies the steady states of system (2.1). First, we determine the steady states, and then study their stability. 4.1. Steady states. Clearly, the origin z0 = (0, 0, 0, 0, 0) in R5 + and zλ = (0, 0, λ, 0, 0) are steady states. As we show below, the origin is unstable while zλ is stable. Then, we turn to find the other possible steady states z = (xc, xe, yr, yn, w), and assume that yr > 0, since without resources the model does not make sense, and x̄c > 0, since without workers there is no wealth generation. To find the steady states, we solve the nonlinear algebraic system 0 = (βc − αc)xc + γexe − γcxc, 0 = (βe − αe)xe − γexe + γcxc, 0 = γyr(λ− yr)− δrxcyr, 0 = yn 0 = δwxcyr − Cc − Ce, (4.1) where Cc and Ce are given in (2.9), wth is given in (2.7), and αc and αe are given in (2.10). Next, we omit the ε in wth, since x̄c > 0, and use wth = ρ(xc + κ0xe). We introduce two nonnegative parameters that simplify the writing in what follows: η = w wth , ψ = xe xc , (4.2) which measure the ratio of wealth to the wealth threshold, and the ratio of Elites to Commoners. Then, Cc = min(1, η)sxc, Ce = min(1, η)κ1sxe, and αc = αm + max(0, 1−min(1, η))(αM − αm), αe = αm + max(0, 1− κ1 min(1, η))(αM − αm). We recall that, by assumption, 0 < αm ≤ βe ≤ βc ≤ αM < 1. We note that κ1 ≥ 1 and consider three cases: (i) η < 1/κ1, in which there is not enough wealth for the minimal existence of both populations; (ii) 1/κ1 ≤ η < 1, in which there is enough wealth for the rich but not enough for the workers; and (iii) 1 ≤ η in which there is enough wealth for both populations. Case (i): η < 1/κ1 – not enough wealth. Then, Cc = ηsxc, Ce = ηκ1sxe, and αc = αM − η(αM − αm), αe = αM − κ1η(αM − αm). Since xc > 0 by assumption, the first equation in (4.1) implies, 0 = βc − αM + η(αM − αm) + γeψ − γc. From the second equation, we get either xe = 0, i.e., ψ = 0 (no Elites), or 0 = βe − αM + κ1η(αM − αm)− γeψ + γc. Adding the two expressions and using simple algebra yields (ψ > 0), η = (2αM − βc − βe) (αM − αm)(κ1 + 1) . (4.3) EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 9 Now, the condition η < 1/κ1 reads 2αM − βc − βe αM − αm < 1 + 1 κ1 . (4.4) It follows that when this inequality does not hold, Case (i) cannot be realized as a nonnegative steady solution. Furthermore, we find that the ratio ψ is determined by the coefficients ψ = 1 γe(κ1 + 1) ((κ1 − 1)αM + (κ1 + 1)γc + βe − κ1βc) . (4.5) We note in passing that ψ = 0, i.e., no Elites, can happen only when γc = 0, so no upward mobility, and βe = βc = αM . We describe this case below. The equation for yr leads to 0 = (γ(λ− yr)− δrxc)yr, and since we assume that yr > 0, we find 0 = γ(λ− yr)− δrxc. Hence, xc = γ δr (λ− yr). The equation for w yields 0 = δwxcyr − ηsxc − ηκ1sxe = (δwyr − ηs− ηκ1sψ)xc, and since xc > 0, 0 = δwyr − ηs− ηκ1sψ. Therefore, yr = ηs δw ψ1. This implies that xc = γ δr (λ− yr) = γ δr ( λ− ηs δw ψ1 ) . Next, xe = ψxc. Finally, w = ηρ(xc + κ0xe), thus w = ηρ(1 + κ0ψ)xc. Moreover, since xc > 0, we find that w > 0 and 0 < γ δr ( λ− ηs δw ψ1 ) , which provides the following bound on ψ, 0 ≤ ψ < 1 κ1 (λδw sη − 1 ) . (4.6) Case (ii): 1/κ1 ≤ η < 1 not enough wealth for the workers, but enough for the rich. Proceeding as above, we find Cc = ηsxc, Ce = ηκ1sxe, αc = αM − η(αM − αm), αe = αm. Similarly, since xc > 0, we obtain from the first equation in (2.1), 0 = βc − αM + η(αM − αm) + γeψ − γc. 10 M. SHILLOR, T. A. KADHIM EJDE-2023/59 The second equation yields either xe = 0, which is the case with ψ = 0, or 0 = (βe − αm − γe)ψ + γc. Thus, ψ is also fixed in this case, ψ = γc αm − βe + γe . (4.7) Therefore, 0 ≤ ψ only when 0 < αm − βe + γe. Otherwise, Case (ii) does not yield any steady states. Then, η is determined as, η = αM − βc αM − αm − (αm − βe)γc (αM − αm)(αm − βe + γe) . (4.8) In this case, the right-hand side must satisfy 1/κ1 ≤ η < 1, otherwise Case (ii) cannot have steady solutions. Proceeding as in Case (i), the equations for yr and w lead to the same expressions but with the new η and ψ. We summarize the findings in Cases (i) and (ii) as follows. Proposition 4.1. Assume that (i) 0 ≤ η < 1/κ1; or (ii) 1/κ1 < η < 1 and 0 < αm − βe + γe. Then, the unique steady state solution z = (xc, xe, yr, yn, w), in which xc > 0, is given by xc = γ δr ( λ− ηs δw ψ1 ) , xe = ψxc, yr = ηs δw ψ1, yn = 0, w = ηρ(1 + κ0ψ)xc. (4.9) In Case (i) η is given in (4.3) and ψ in (4.6); while in Case (ii) η is given in (4.8) and ψ in (4.7). Case (iii): 1 ≤ η there is enough wealth for everyone. We note that in this case the equations for xc and xe are not coupled to those for y and w. Proceeding as above, we obtain Cc = sxc, Ce = κ1sxe, αc = αe = αm. Since xc > 0, the first equation in (2.1) leads to 0 = βc − αm + γeψ − γc. From the second equation, we get either xe = 0, which is the case with ψ = 0, or 0 = (βe − αm − γe)ψ + γc. Hence, ψ = γc αm − βe + γe = 1 γe (αm − βc + γc). (4.10) Therefore, 0 ≤ ψ only when 0 < αm − βe + γe and 0 < αm − βc + γc. Otherwise, Case (iii) doesn’t yield any steady states. Then, the equation for yr yields 0 = γ(λ− yr)− δrxc. Thus, xc = γ δr (λ− yr). EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 11 Finally, the equation for w yields 0 = δwyr − s− κ1sψ. Therefore, yr = s δw ψ1. This implies that xc = γ δr (λ− yr) = γ δr ( λ− s δw ψ1 ) . Next, xe = ψxc. Finally, w = ηρ(xc + κ0xe), thus w = ηρ(1 + κ0ψ)xc. Now, by assumption xc > 0, hence w > 0 and 0 < γ δr ( λ− s δw ψ1 ) , which provides the following bound on ψ, 0 ≤ ψ < 1 κ1 (λδw s − 1 ) . We summarize Case (iii) as follows. Proposition 4.2. Assume that 1 ≤ η and 0 < αm − βe + γe, 0 < αm − βc + γc. Then, the unique steady state solution z = (xc, xe, yr, yn, w), is given by xc = γ δr ( λ− s δw ψ1 ) , xe = ψxc, yr = s δw ψ1, yn = 0, w = ηρ(1 + κ0ψ)xc. (4.11) Here, ψ is given in (4.7) and η in (4.8). Remark 4.3. We note that when 1 ≤ η only the wealth w depends on η, while ψ is fixed by the system parameters. Thus, (xc, xe, yr, yn) are uniquely determined, so they have the same values, independently of the wealth, as long as 1 ≤ η, i.e., there is sufficient wealth to have the populations be above the famine threshold. 4.2. Stability. We turn to study the stability of the steady states. To that end, we use the Jacobian matrix J of F ( constructed in the Appendix) evaluated at the steady states. It is found that we have five versions of J corresponding to the two ‘simple’ steady states z = (0, 0, 0, 0, 0) and z = (0, 0, λ, 0, 0), and three different cases of the steady states derived above. To simplify the expressions, we use the notation ∆α = αM − αm, ψ0 = (1 + κ0ψ)−1, ψ1 = 1 + κ1ψ. In the ‘simple’ steady states, we have αc = αe = αM , hence ∆α = 0; η = 0 since w = 0; and ψ = 1, hence ψ0 = 1/(1 + κ0) and ψ1 = 1 + κ1, then it is found 12 M. SHILLOR, T. A. KADHIM EJDE-2023/59 (Appendix) that when we substitute z = (0, 0, 0, 0, 0) into J(i), we obtain J(0, 0, 0, 0, 0) =  βc − αM − γc γe 0 0 0 γc βe − αM − γe 0 0 0 0 0 γλ 0 0 0 0 0 − δnγλδr 0 s(1 + κ1) 0 δwγλ δr δwγλ δr − s(1+κ1) ρ(1+κ0)  . (4.12) The eigenvalues are λ1,2 = 1 2 ( βc − 2αM − γc + βe − γe ± √ ∆ ) , λ3 = γλ, λ4 = −δnγλ δr , λ5 = − s(1 + κ1) ρ(1 + κ0) , where ∆ = (βc − αM − γc)2 − 2(βc − αM − γc)(βe − αM − γe) + (βe − αM − γe)2 + 4γcγe. We note that λ3 is positive, hence the origin is unstable, actually, repelling. Next, a similar calculation shows that J(0, 0, λ, 0, 0) =  βc − αM − γc γe 0 0 0 γc βe − αM − γe 0 0 0 0 0 −γλ 0 0 0 0 0 − δnγλδr 0 s(1 + κ1) 0 δwγλ δr δwγλ δr − s(1+κ1) ρ(1+κ0)  . (4.13) Then, the eigenvalues λ1,2 are as above and because of the assumption that 0 ≤ ψ only when 0 < αm − βe + γe and 0 < αm − βc + γc, they have negative real parts. Also, λ4 and λ5 are the same as above while λ3 = −γλ < 0. Therefore, the steady state is stable and attracting. We note that in the original HANDY model, [5], κ0 = κ1 = κ, δr = δw = δ, δn = 0, and since there is no social mobility, γc = γe = 0. It follows that in their case, too, the origin is unstable while the steady state (0, 0, λ, 0) is stable and attracting. These results show that the model allows for the collapse of the society, so that only nature remains, while the populations decay to zero. We turn to the other steady states. To simplify the presentation, we use the notation ∆α = αM − αm, ψ0 = 1 1 + κ0ψ , ψ1 = 1 + κ1ψ. In Case (i) we find J(i)(z) =  (βc − αc − γc) − η∆αψ0 γe − ηκ0∆αψ0 0 0 ∆αψ0ρ −1 γc − ηκ1∆αψψ0 J22 0 0 κ1∆αψψ0ρ −1 − ηsδr δw ψ1 0 J33 0 0 0 0 0 −δnγ δr ( λ− s δw ψ1 ) 0 sψ1 − ηs(κ0 − κ1)ψ0ψ ηs(κ0 − κ1)ψ0 J53 γδw δr ( λ− ηs δw ψ1 ) −sψ1ψ0ρ −1  (4.14) EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 13 Where, J22 = (βe−αe−γe)−ηκ0κ1∆αψ0ψ, J33 = −γηsδw ψ1, and J53 = γδw δr ( λ− ηs δw ψ1 ) . Similarly, in Case (ii), we have J(ii)(z) =  (βc − αc − γc) − η∆αψ0 γe − ηκ0∆αψ0 0 0 ∆αψ0ρ −1 γc βe − αe − γe 0 0 0 − ηsδr δw ψ1 0 0 − γηs δw ψ1 0 0 0 0 −δnγ δr ( λ− s δw ψ1 ) 0 sψ1 − ηs(κ0 − κ1)ψ0ψ ηs(κ0 − κ1)ψ0 J53 γδw δr ( λ− ηs δw ψ1 ) −sψ0ψ1ρ −1  (4.15) Where, J53 = γδw δr ( λ− ηs δw ψ1 ) . Finally, in Case (iii), we obtain J(iii)(z̄) =  βc − αm − γc γe 0 0 0 γc βe − αm − γe 0 0 0 − sδrδw ψ1 0 − sγ δw ψ1 0 0 0 0 0 −δnγ δr ( λ− s δw ψ1 ) 0 κ1sψ −κ1s γδw δr ( λ− s δw ψ1 ) γδw δr ( λ− s δw ψ1 ) 0  . (4.16) Clearly, the study of the states stability can be done only in particular cases when the constants are given. Thus, we have closed form expressions for the Jacobian matrices in the five cases. In each of the simulations below, we also included the information about which case is simulated and its stability, when it converges to a steady state. We note that in the baseline simulations of Case (iii), with a modified death rate so that η ≥ 1, we found λ1 = 0, λ2 = −0.76 λ3 = −1.67, λ4 = 0.02, λ5 = −0.07. This means that the steady state is unstable, and further study may show periodic solutions, non-periodic oscillations or chaotic oscillations. 5. Simulations This section presents and briefly discusses the results of our model numerical experiments. The main interest is to exhibit the four different types of behavior of the solutions: monotone or oscillatory approach to the steady state; periodic solutions; and system collapse. It is very likely that the model can predicts chaotic behavior, however, we did not observe chaos with the system parameters we have experimented with. The description of the code for the simulations can be found in [1]. We start with a short study of the effects of increasing δn, the the renewables depletion coefficient, and the effects of doubling κ1, the pay gap measure. Table 1 lists the model dependent variables, their initial values, and the values of the parameters. These were chosen in the baseline simulation, which depicts a typical system behavior. To obtain the four different types of behavior, we changed some of the parameters and initial conditions, and the changes are described with each type of behavior. The results are depicted in the Figs. 1–6. The numbers in the table include multiples of the renewable resources depletion coefficient, to highlight the over consumption of resources. 14 M. SHILLOR, T. A. KADHIM EJDE-2023/59 In all the simulations, the blue curve depicts the workers xc; magenta the rich xe; green and black curves are the renewables yr and nonrenewables yn, respectively; and the society’s wealth w is represented by the cyan curve. We note that here, unlike [5], we use ‘time units’ instead of ‘years,’ on the horizontal scale, since our units are quite arbitrary, so the time is not necessarily identifiable as a day, month or year. The main purpose is to depict the different types of solutions. Table 1. Parameters and data used in the simulations of the HANDY-SMRN model. Symbols Values Meaning ∆t 0.05 time unit τF 1000 and 10,000 time period xc(0) 10 Commoners initial condition xe(0) 2 Elites initial condition yr(0) 100 initial Renewables yn(0) 100 initial Non-Renewables w(0) 0 initial wealth λ 100 Natures carrying capacity γ 0.01 Natures regeneration factor δ∗ or δr 6.66 · 10−6 and multiples Natures depletion factor δn 6.15 · 10−6 Nonrenewables depletion factor αM 0.07 Maximum expected death rate αm 0.01 Normal death rate βc 0.03 Commoners birth rate βe 0.03 Elites birth rate s 5 · 10−4 salary per capita ρ 1 · 10−4 minimum required consumption κ0 1 wealth unbalanced threshold rate κ1 1 or 2 Elites rate of consumption γc 0.085 or 0.075 Commoners to Elites mobility factor γe 0.075 or 0.085 Elites to Commoners mobility factor δw 1.18 · δr wealth rate of growth 5.1. Qualitative behavior. The four model qualitatively different types of so- lutions are depicted in Figure 1, and all fit Case (i) of the analysis, in which η ≤ 1/κ1 = 1, since we use κ1 = 1. Figure 1 (TL-top left) describes a ‘soft landing’ [7] in which the system ap- proaches the steady state of (45000, 4000, 60, 0, 5.5) monotonically after a single initial oscillation. From time 500 onward, the system converges monotonically to this steady state. The nonrenewables tend to zero, as expected, while the renew- able resources tend to 60, which is lower than λ = 100. The two populations and the wealth converge to their limits smoothly. The depletion factor δr = 1.3 δ∗ has a different value than in Table 1. Also, η = 0.67 at the steady state. The two populations grow to high values causing complete depletion of the nonrenewables and almost all of the wealth, thus, leaving the society’s populations with only the EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 15 Figure 1. Four basic scenarios of long-time behavior and ap- proach to a steady state in the HANDY-SMRN model: (TL) mono- tone approach after initial oscillation; (TR) decreasing oscillations; (BL) steady oscillations; (BR) collapse; κ1 = 1, η < 1 = 1/κ1 (Case (i) above). renewables. The wealth decrease to 5.5 from its full value of 20 at the peek of the population rise. The vanishing of the nonrenewables and a considerable decrease in wealth drag the populations to half of their prosperity peeks. Figure 1 (TR-top right) depicts convergence to the steady state after a long period of large but damped oscillations, and then after time 575, the oscillations al- most vanish and the steady state is rapidly approached. We note that, as expected, the nonrenewables do not oscillate, since these are just consumed. The value of δr = 4.6 δ∗ was different than that in Table 1, and we increased the time span for a better depiction of the approach to the steady state. The steady state the system settles to around (27000, 25000, 20, 0, 0). Moreover, η = 0.67 at the steady state, which is Case (i). Figure 1 (BL bottom left) provides numerical evidence that the system may have periodic solutions in (xc, xe, yr, w), once the initial conditions have been cho- sen appropriately. However, because of the structure of the system, the nonre- newables yn cannot oscillate, and must decay. Here, δr = 7.2δ∗, and a shift of the initial conditions was made to (750, 750, 96, 0, 0) to generate the periodic solution. Moreover; the following changes were also used, (γec , γ c e , δr, δn, δw) = (0.085, 0.075, 7.2 δ∗, 6.15 · 10−6, 1.18δr). The mathematical proof of the existence of periodic solutions is unavailable, yet, however, these simulations provide a strong indication that they are likely to exist. Finally, Figure 1 (BR) depicts the system collapse at about time 600. The populations increase considerably from the initial conditions and then the resources, the wealth and the populations collapse. We note that the populations’ values at nature’s collapse are of interest. For instance; in Figure 1 (BL) the renewables and wealth collapse around time 210. The populations collapse after another hundred time units, when the commoners population decreases to 850; however, the system starts to recover after another 100 time units, and the renewables increase about time 310 with the commoners 16 M. SHILLOR, T. A. KADHIM EJDE-2023/59 working to revive the wealth. In the meantime; when δr is large, such as in Figure 1 (BR), the excessive consumption of resources doesn’t allow the system to recover; keeping η close to zero which is an unstable case for the system. 5.2. Dependence on δn, δr. We study numerically the system’s dependence on the parameters δn and then δr, varying the rates of depletion and consumption of the nonrenewables and renewables, respectively. This is the new part of the modified model. When the depletion factor of the nonrenewables is decreased, δn = 6.15 · 10−8, i.e. 100-fold less than δn = 6.15 · 10−6, which means a more moderate use of nonrenewables, we find, Figure 2(L), that the system reaches a steady state after 2000 time units, and the populations’ numbers are considerably larger than in the case of increased depletion, δn = 6.15 · 10−5, Figure 2 (R). In both cases the system eventually reaches its steady state, but, in both cases the populations almost collapse following the collapse of the renewables. We mention that a decrease in δn causes the wealth to increase; that is also the reason behind the expansion of the two populations. That is not the case with an increased δn, shown in Figure 2 (R), since the nonrenewables are drained because of the overuse in a shorter time period, so the system depends only on the renewables as a source of wealth. In both figures the change was done only to δn; while keeping the constants at the values (γec , γ c e , δr, κ1, δw) = (0.085, 0.075, 2.19 δ∗, 1, 1.18δr). Figure 2. Two scenarios with different δn; (L) decreased δn; (R) increased δn; κ1 = 1 (Case (i) of the analysis above). To obtain some insight into the effects of increasing the use of renewable re- sources, we performed an experiment by changing the renewables rate constant δr, and the results are provided in Figures 3–5. Here, we also use κ1 = 2. Figure 3. A scenario with δr = 9.9δ∗; (L) Fast damped oscilla- tions; (R) detail zoom; κ1 = 2 (Case (i) in the analysis). It is found that within the whole interval 9.96δ∗ ≤ δr ≤ 23.2δ∗ there are rapid oscillations. However, Figure 3 (L), where δr = 9.9δ∗, seems to be of interest, since there are many fast slightly damped oscillations until about time 7100, and then, EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 17 in a very short period of about 100 time units, the system approaches the steady state. The solution about time 7000, and the fast approach to the steady state, are depicted in the zoom 3 (R). On the other hand, in Figure 4 we use δr = 23.2δ∗ and it seems that the solutions are periodic. We use the same parameters in both Figures 3–5, except for δr. Figure 4. A scenario with δr = 23.2δ∗; (L) the system settles quickly into periodic oscillations with period of about 1000 time units; (R) zoom detail. Figure 5. The scenario with δr = 23.2δ∗, (L) consumption rates and (R) death rates. For the sake of completeness, in a shorter time period of 5000 time, Figure 5 depicts the death and consumption rates as functions of time for the scenario with δr = 23.2δ∗. Moreover, it is found that during the oscillations, the system approaches partial collapse and then it recovers. We note an interesting behavior in Figure 5 (L), since having piecewise constant values in a dynamical system, for the consumption rates, is somewhat unusual, and here it seems to be related to the non-differentiability of the rate functions. Also, the death rates have rather interesting spikes. 5.3. Dependence on κ1. Finally, we study the system’s dependence on the con- sumption ratio κ1; noting that in the previous simulations, we used κ1 = 1. Here, we used four cases κ1 = 0.5, 1, 5, 10. The simulations results are depicted in Figure 6. We use γe = 0.075 and γc = 0.085. The case κ1 = 0.5 is both unrealistic in terms of the model, since the rich would not consume less than the commoners, and mathematically it produces what seems to be a periodic solution. The other three cases lead to very similar approach to the steady states and comparable steady states, (z = xc, xe, yr, yn, w). The convergence is either monotone or with a few 18 M. SHILLOR, T. A. KADHIM EJDE-2023/59 oscillations before the monotone convergence. Indeed, the steady states were found to be (approximately), κ1 = 1, z = (29000, 33000, 23, 0, 4); κ1 = 5, z = (27500, 41600, 27, 0, 1); κ1 = 10, z = (22600, 35000, 40, 1). However, we note that the ways they were reached seem to be quite different, with more or less societal upheaval in the form of large oscillations. This indicates that just studying the steady states of a society’s model may not provide a good description and insight into the society’s dynamics. Figure 6. Four scenarios with κ1 = 0.5 (TL), κ1 = 1(TR), κ1 = 5 (BL), κ1 = 10 (BR), in all cases δr = 4δ∗ 6. Conclusions and further research This work presents the HANDY-SMRN model, which is a new version of the HANDY model [5], that is more realistic. It is a continuation of the studies of other versions of the HANDY model in [1]. This HANDY-SMRN model extends the HANDY model by taking into account mobility between the Commoners (work- ers) and Elites (rich), and splits the natural resources into renewables and nonre- newables, where the latter once consumed cannot be renewed. The model is a nonlinear coupled system of five ODEs for the Commoners and Elite populations, societal wealth, and renewable and nonrenewable resources. It establishes the ex- istence, positivity and boundedness of the solutions of the model. We note that similar properties of the HANDY model have been established in [1], and the anal- ysis here is similar, but tailored to the HANDY-SMRN model. Then, it studies the system’s steady states and their stability. It is found that the steady states depend essentially on the ratio η = w/wth of the wealth w to the wealth’s threshold wth. This type of models, which describe in a very coarse way complex societies, although relatively ‘simple,’ seem to provide insight into what may be common in the collapse of such societies. Clearly, there is considerable interest to find out if the model can fit some of the historical collapses of societies which essentially consisted of commoners and nobility. EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 19 Mathematically, there is interest in establishing that the model has oscillatory solutions (except for the nonrenewables), since our computer experiments indicate that this is the case. Moreover, it is very likely that a fifth type of behavior exists, namely chaotic behavior, and it is of interest to obtain numerically chaotic solutions, and possibly to prove the existence of chaos rigorously. However, these questions remain unresolved, yet, and will need future investigations. The area is wide open and there is plenty of room for various types of mathe- matics that may shed light on social dynamics. 7. Appendix We compute the Jacobian J of F, proceeding entry by entry, then evaluate it at the corresponding steady states z0 = (0, 0, 0, 0, 0), zλ = (0, 0, λ, 0, 0), and Cases (i)–(iii). J11 = ∂F1 ∂xc = (βc − αc − γc)− xc ∂αc ∂xc , J12 = ∂F1 ∂xe = γe − xc ∂αc ∂xe , J13 = ∂F1 ∂yr = 0, J14 = ∂F1 ∂yn = 0, J15 = ∂F1 ∂w = −xc ∂αc ∂w . J21 = ∂F2 ∂xc = γc − xe ∂αe ∂xc , J22 = ∂F2 ∂xe = (βe − αe − γe)− xe ∂αe ∂xe , J23 = ∂F2 ∂yr = 0, J24 = ∂F2 ∂yn = 0, J25 = ∂F2 ∂w = −xe ∂αe ∂w . J31 = ∂F3 ∂xc = −δryr; J32 = ∂F3 ∂xe = 0, J33 = ∂F3 ∂yr = γλ− 2γyr − δrxc, J34 = ∂F3 ∂yn = 0 J35 = ∂F3 ∂w = 0. J41 = ∂F4 ∂xc = −δnyn, J42 = ∂F4 ∂xe = 0, J43 = ∂F4 ∂yr = 0, J44 = ∂F4 ∂yn = −δnxc, J45 = ∂F4 ∂w = 0. Finally, J51 = ∂F5 ∂xc = δw(yr + yn)− ∂Cc ∂xc − ∂Ce ∂xc , J52 = ∂F5 ∂xe = −∂Cc ∂xe − ∂Ce ∂xe , J53 = ∂F5 ∂yr = δwxc, J54 = ∂F5 ∂yn = δwxc, J55 = ∂F5 ∂w = −∂Cc ∂w − ∂Ce ∂w . To compute the partial derivatives of αc, αe, Cc, Ce, we consider the five cases. To simplify the notations below, we use the equilibrium solutions, and let ∆α = αM − αm, ψ0 = (1 + κ0ψ)−1, ψ1 = (1 + κ1ψ), and recalling that η = w/ρ(xc + κ0xe + ε) (where we use ε = 0, unless we need to guarantee that the denominator does not vanish). We begin with the Jacobian evaluated at z0 = (0, 0, 0, 0, 0), zλ = (0, 0, λ, 0, 0) and Case (i). To that end we note that in all three cases 0 ≤ η < 1/κ1, and so, Cc = ηsxc, Ce = ηκ1sxe, αc = αM −η(αM −αm), αe = αM −κ1η(αM −αm). 20 M. SHILLOR, T. A. KADHIM EJDE-2023/59 Therefore, ∂αc ∂xc = −∆α ∂η ∂xc = η∆αψ0 1 xc , ∂αc ∂xe = ηκ0∆αψ0 1 xc , ∂αe ∂xc = −κ1∆α ∂η ∂xc = ηκ1∆αψ0 1 xc , ∂αe ∂xe = ηκ0κ1∆αψ0 1 xc , ∂αc ∂w = −∆αψ0 1 ρxc , ∂αe ∂w = −κ1∆αψ0 1 ρxc , ∂Cc ∂xc = ηs+ sxc ∂η ∂xc = ηsκ0ψ0ψ, ∂Cc ∂xe = sxc ∂η ∂xe = −ηsκ0ψ0, ∂Cc ∂w = sψ0 ρ , ∂Ce ∂w = sκ1ψψ0 1 ρ , ∂Ce ∂xc = κ1sxe ∂η ∂xc = −ηsκ1ψψ0; ∂Ce ∂xe = ηsκ1ψ0. Then, at the solution z0 = (0, 0, 0, 0, 0), we have, η = 0 (since w = 0 and we use ε > 0), and then J11 = βc − αM − γc, J12 = γe, J13 = J14 = J15 = 0. J21 = γc, J22 = βe − αM − γe, J23 = J24 = J25 = 0. J31 = J32 = 0, J33 = γλ, J34 = J35 = 0. J41 = J42 = J43 = 0, J44 = −δnγλ δr , J45 = 0. J51 = s(1 + κ1), J52 = 0, J53 = J54 = γδwλ δr , J55 = − s(1 + κ1) ρ(1 + κ0) . The matrix is given now in (4.12). We perform a similar computation for z0 = (0, 0, λ, 0, 0), and the resulting matrix is given in (4.13). We next consider Case (i) when 0 ≤ η < 1/κ1. Using these results yields the following Jacobian, evaluated at a steady state of Case (i), (Proposition 4.1). J11 = (βc − αc − γc)− η∆αψ0, J12 = γe − ηκ0∆αψ0, J13 = 0, J14 = 0, J15 = ∆αψ0ρ −1. Next, J21 = γc − ηκ1∆αψψ0, J22 = (βe − αe − γe)− ηκ0κ1ψ0∆αψ, J23 = 0, J24 = 0, J25 = κ1∆αψψ0ρ −1. Furthermore, J31 = −δryr = −ηsδr δw ψ1, J32 = 0, J33 = −γηs δw ψ1, J34 = 0, J35 = 0. J41 = −δnyn = 0, J42 = 0, J43 = 0, J44 = −δnγ δr ( λ− ηs δw ψ1 ) , J45 = 0. Finally, J51 = sψ1 + ηs(κ1 − κ0)ψψ0, J52 = ηs(κ0 − κ1)ψ0, J53 = J54 = γδw δr ( λ− ηs δw ψ1 ) , J55 = − s ρ ψ1ψ0. EJDE-2023/59 ANALYSIS AND SIMULATIONS OF THE HANDY MODEL 21 The Jacobian matrix in Case (i), with the entries above evaluated at the steady state (Proposition 4.1), is given by (4.14). We consider now Case (ii), where Cc = ηsxc, Ce = ηκ1sxe, αc = αM − η(αM − αm), αe = αm. It follows, since αe = αm, that the only changes from Case (i) are those related to αe. Thus, we find that the only differences in (4.14) are in the following components: J21 = γc, J22 = βe − αm − γe, J25 = 0, and this yields (4.15). Finally, in Case (iii), we have Cc = sxc, Ce = κ1sxe, αc = αe = αm. Then, the changes in J(z̄), as compared to Case (i) are those related to αc and αe. Therefore, ∂Cc ∂xc = s, ∂Cc ∂xe = 0, ∂Ce ∂xc = 0, ∂Ce ∂xe = κ1s, ∂Cc ∂w = ∂Ce ∂w = 0. Thus, the changes in the entries, using the steady values given in Proposition 4.2, are the following: J11 = βc − αm − γc, J12 = γe, J13 = 0, J14 = 0, J15 = 0. J21 = γc, J22 = βe − αm − γe, J23 = 0, J24 = 0, J25 = 0. J31 = −sδr δw ψ1, J32 = 0, J33 = −γs δw ψ1, J34 = 0, J35 = 0. J41 = 0, J42 = 0, J43 = 0, J44 = −δnγ δr ( λ− s δw ψ1 ) , J45 = 0. J51 = δw yr − s = κ1sψ, J52 = −κ1s, J53 = J54 = γδw δy ( λ− s δw ψ1 ) , J55 = 0. Assembling J(iii) in Case (iii) yields (4.16). Acknowledgments. The authors would like to thank the reviewer for useful com- ments. References [1] T. A. K. A. Al-Khawaja; Mathematical Models, Analysis and Simulations of the HANDY Model with Middle Class, Doctoral dissertation, Oakland University, 2021. [2] J. A. Brander, M. S. Taylor; The simple economics of Easter Island: A Ricardo-Malthus model of renewable resource use, American economic review, textbf88(10 (1998), 119–138. [3] B. Grammaticos, R. Willox, J. Satsuma; Revisiting the Human and Nature Dynamics Model. Regular and Chaotic Dynamics, 25(2) (2020), 178–198. [4] K. Kuttler; Modern Analysis, CRC press, 1997. [5] S. Motesharrei; Minimal Models of Human-Nature Interaction, Doctoral dissertation, Uni- versity of Maryland, 2014. [6] S. Motesharrei, J. Rivas, E. Kalnay; Human and nature dynamics (HANDY): Modeling inequality and use of resources in the collapse or sustainability of societies, Ecological Eco- nomics, 101 (2014), 90–102. [7] S. Motesharrei, J. Rivas, E. Kalnay, et al.; Modeling sustainability: population, inequality, consumption, and bidirectional coupling of the Earth and Human Systems, National Science Review, 3(4) (2016), 470–494. Meir Shillor Department of Mathematics and Statistics, Oakland University, Rochester, MI 48309, USA Email address: shillor@oakland.edu 22 M. SHILLOR, T. A. KADHIM EJDE-2023/59 Thanaa Ali Kadhim Department of Mathematics and Statistics, Oakland University, Rochester, MI 48309, USA Email address: taalkhawaja@oakland.edu 1. Introduction 2. HANDY-SMRN model 3. Model analysis 3.1. Positivity and a priori estimates 4. Steady states and stability 4.1. Steady states 4.2. Stability 5. Simulations 5.1. Qualitative behavior 5.2. Dependence on n, r 5.3. Dependence on 1 6. Conclusions and further research 7. Appendix Acknowledgments References