A Nonlinear ODE Model for a Consumeristic Society

: In this paper, we introduce an ODE system to model the interaction between natural resources and human exploitation in a rich consumeristic society. In this model, the rate of change in population depends on wealth per capita, and the rate of consumption has a quadratic growth with respect to population and wealth. We distinguish between renewable and non-renewable resources, and we introduce a replenishment term for non-renewable resources. We first obtain some information on the asymptotic behavior of wealth and population, then we compute all system equilibria and study their stability when the resource exploitation parameter is low. Some numerical simulations confirm the theoretical results and suggest directions for future research.


Introduction
In this paper, we study a model of interactions between society and nature, whose original inspiration goes back to the HANDY model introduced in [1] (see also [2,3]).The original HANDY model is a set of four differential equations whose variables are natural resources, wealth produced by human work, and population, split in the two classes of Commoner and Elite.Using computer simulations, the paper [1] shows different possibilities for the evolution of society, including collapse.It is a simple model that seems to allow a rich panel of behavior, so it attracted attention, and several papers have been devoted to its exploration, trying to generalize the original idea in several, different directions.In [4], the HANDY model is studied using a Lyapunov-type analysis.The authors give an interesting mathematical generalization of the original model, proving that a set of qualitative hypotheses on the state variables implies collapse, even without assuming that the functions involved are solutions of a system of differential equations.In [5], the authors on one hand simplify the HANDY model, dropping the division of population in two classes, and on the other hand use general nonlinear functions to describe the equations involving the state variables, which are "population", "resources", and "reserves".Because of the generality of their model, they are able to treat different sets of hypotheses on the interactions between the state variables.More recently, Tonnelier [6] presented a clearer rewrite of the HANDY model and studied its dynamics using bifurcation techniques and focusing on the influence of two parameters: the nature depletion rate and the inequality factor.All these works, after the mathematical analysis, give several kinds of numerical simulations, exploring the actual dynamics of the system for some fixed values of the parameters.Computational techniques applied to the HANDY model are the main subject of [7].
Our work is more directly linked to some ideas introduced in [8,9].In [8], various developments of the original model are introduced, including, among other aspects, the division of natural resources into renewable and non-renewable categories.This idea has been developed in the paper [9], where some general results are obtained and again numerical simulations show different possibilities for the evolution of society.In [9] the dynamics of non-renewable resources is given by the following equation where y n is the variable for the non-renewable resources, x c is the Commoner population and δ n > 0 is a parameter.Clearly, this equation translates the idea of a given amount and an irreversible depletion of y n .In our previous work [2], we have pursued this idea but we have introduced a replenishment term for non-renewable resources.This is due to the optimistic argument, which is often used in public debate on these subjects, that human ingenuity and scientific progress can substitute depleted resources with new ones.The term of replenishment that we have added is given by This term depends on x w, where x represents the population, w represents wealth, z is the level of non-renewable resources, and k is a constant.Hence, the dynamics of non-renewable resources z in [2] is given by The replenishment term is obtained from a function t t+1 which, in population dynamics, is a Holling II-type function (see [10], p. 25 and passim), used to model saturation effects.The equation above states that the replenishment is possible but cannot be above a level k, and we can call this an hypothesis of "moderate optimism".Another characteristic of [2] is that we dropped the distinction between Commoner and Elite, so we obtained a model with four variables, which are population, renewable and non-renewable resources, wealth.In [2], we obtained all the critical points of the system, and we studied their stability.In the final section of [2], numerical simulations support the theoretical results and also suggest the possibility, not yet treated in a rigorous way, of periodic orbits and chaotic dynamics.
The present paper originated from the attempt to adapt the HANDY model to contemporary, rich, consumeristic societies (like contemporary Western countries).Indeed, the HANDY model seems to be a very nice and simple way to treat the interaction of nature and society for a large class of historical societies, but, in our opinion, it has some serious flaws when dealing with contemporary, rich, consumeristic societies.The main one is probably the fact that in the original model, the birth rate is constant and the death rate decays when wealth increases, leading to the conclusion that a wealthy society should have a strong population growth.This is definitely not true in contemporary Western societies, and the reason for this is due to the fact that birth rate is not constant, and indeed, in most contemporary rich societies, it is falling.Hence, to study rich, contemporary societies, it seems necessary to overcome the HANDY model used in the papers quoted above (including our previous paper [2]) and to introduce new ideas.This is what we do in the present paper.Firstly, we introduce an equation for population dynamics which is different from that of HANDY model, and it is as follows Here, x is the population and w the wealth.In this equation we have a negative term −cx, which is balanced by the term depending on the wealth per capita w/x via a constant µ.When w/µx is very small or very large, the positive term is small and the population falls.In an intermediate region, the population grows.This seems to better fit the recent history of rich countries.
Another difference compared to the HANDY model (and to [2]) is that in previous works the depletion of natural resources depends on the population, while it seems reasonable that in a consumeristic society it should depend more on wealth than on population.Hence, a first idea should be to substitute population with wealth (x with w) in the depletion terms, giving rise to quadratic terms w 2 .We will pursue this idea in a forthcoming paper.In the present paper, however, we are interested in preserving some influence of population on consumption of natural resources, so we introduce the function m = m(x, w) = max{x, w} and we substitute x with m in the depletion terms.We notice also that the consumption term in the equation for w ′ in [2], which is given by s x − 1 ρ (ρx − w) + , has a linear growth, while here, to model a consumeristic society, we use a quadratic term m(x, w)w, where the consumption is led mostly by wealth.
Summing up all these ideas, and applying some standard rescaling, the model we are going to study in the present paper is the following where x is the population, y and z are the renewable and non-renewable resources, w is the wealth.As we said above, m(x, w) = max{x, w} = 1 2 (x + w + |x − w|).We underline that, because of the changes that have been explained above, the present model can no longer be considered a HANDY-type model, as it was the model of our previous work [2].
To start the work with (1) we need some assumptions on the parameters.Firstly, we assume d > 2c, which is necessary to have a positive growth of x for some range of w µx .
for any value of x, w, hence it is always x ′ ≤ 0, making the model less interesting.We also assume µ > 1, which means that the ratio w/µx becomes large only when w is larger than x, and this implies that the negative effect on the population growth of the increasing wealth becomes relevant only for large values of wealth per capita.We also assume µ ̸ = d+ , which is just a technical assumption that we will drop in a forthcoming paper.We remark that in (1), the depletion rates for renewable and non-renewable resources are equal.This is a simplifying assumption that we hope to drop in future research.In our previous paper [2], numerical simulations seemed to suggest that the dynamics do not change too much taking different values for the depletion rates.In Section 7, some simulations show a different dependence of the dynamics on the two depletion rates (which are labeled δ 1 and δ 2 ), and this seems an interesting topic for future research.
The present paper is organized as follows.After the introduction, in Section 2 we obtain general results for the existence and positivity of the solutions, as well as results on their asymptotic behavior.In particular, we prove that it cannot happen x → +∞ or w → +∞ or z → +∞ as t → +∞.
In Sections 3-6, we compute all the equilibrium points of system (1), and we study their stability properties.We distinguish the cases x < w, x > w, x = w = 0, while the case x = w > 0 is excluded by our hypotheses.There is a rich variety of such equilibria.We are able to give a complete description of the stability properties for some of them, but it seems very difficult to do the same for all.So we decided to study the stability when δ → 0 + , which is reasonable, in our opinion, because δ is a parameter depending on human choices.Using standard linearization techniques and asymptotic developments, we are able to describe the stability properties of all the equilibria when δ → 0 + (and the other parameters are fixed).
In Section 7, we present some numerical simulations, both to corroborate the theoretical results of the preceding sections and to have some hint on the cases for which we do not have such results.
In Section 8 we highlight the main results of the paper, and give some hints on possible future research work.

General Results on the Solutions
We set X = (x, y, z, w) ∈ R 4 and We are interested in non-negative solutions, so we work the cone or in the set Let us first prove existence and uniqueness.
Proposition 1.For any (t 0 , X 0 ) ∈ R × C 0 , there exists a unique maximal solution X(t) to the Cauchy problem X defined in J = (a, b), with a < t 0 < b.
Proof.The function m(x, w) is locally Lipschitz in all R 4 .As a consequence, for any X 0 ∈ C 0 the vector field F(X) is Lipschitz in a neighborhood of X 0 , hence the result follows from standard ODE theory.
In the following, we will set t 0 = 0 and X 0 ∈ C 0 .
The first thing we prove is, if X 0 is as above, the solution stays in C for any t > 0.
Proof.It follows that x(t), y(t), z(t) > 0 for any t ∈ J = (a, b), because the zero constant is a solution in the three equations for x ′ (t), y ′ (t), z ′ (t), so it cannot be crossed, and x 0 , y 0 , z 0 > 0.
Proof.We argue by contradiction, so we assume b < +∞.From the first equation, we easily obtain x ′ ≤ c 1 x for some c 1 > 0, hence As for y(t), using the same argument of [2], we obtain that y is bounded in [0, b) and, more precisely, we have 0 ≤ y(t) ≤ max{λ, y 0 }, ∀t ∈ [0, b).Notice that this also holds when b = +∞.
From the equation for z(t) we obtain z ′ ≤ kz, hence z(t) ≤ z 0 e kt ≤ z 0 e kb := c 3 , and thus, From this it follows that w(t) ≤ c 6 if b < +∞.Therefore, if b < +∞ we obtain that X(t) is bounded in [0, b), which is impossible because of well-known ODE theorems.Hence, b = +∞.
We now prove some asymptotic results on x and w.Proposition 4. It cannot be lim t→+∞ x(t) = +∞.
Proof.We argue by contradiction, so we assume lim t→+∞ x(t) = +∞.For the sake of simplicity, let us set m(t) = m(x(t), w(t)).As x(t) → +∞, of course we have also m(t) → + ∞.To obtain the result, we have to study the asymptotic behavior of y, z, w, in the hypothesis x(t) → +∞.
Proof.Since m = m(x, w) ≥ w we can repeat the previous arguments and obtain lim t→+∞ y(t) = lim t→+∞ z(t) = 0. From this we deduce, again as above, that w(t) is bounded, and this is a contradiction.Proposition 6.It cannot be lim t→+∞ z(t) = +∞.

Steady States If x ≤ w
We are now going to compute all the steady states of the system (1) in C, and to study their stability.As we said in the introduction, in most cases, we will only be able to obtain an answer about the stability for δ → 0 + .In this section we deal with the case x ≤ w.In this subset the equilibrium points are the solutions of the following system w[δy We can consider the following cases: (i) w = 0.In this case we obtain x = 0, and from (4) we have two possibilities: either y = 0 or y = λ.There is no constraint on z, so we have two families of critical points: {(0, λ, Z, 0), Z ≥ 0}.
We remark that the points (7) correspond to what is called "desert state" in [6], while the points (8) correspond to what is called "nature state" in the same paper.We also notice that these points are not in C 0 , and the vector field F is not defined on them.However, it is interesting to study what we can improperly call their stability, that is to know whether the trajectories starting close to them are attracted or repelled by them, because this says to us if a society with declining population and wealth can recover or not.
(ii) w ̸ = 0.In this case we have and we can count 8 different cases.We note that, if x ̸ = 0 from (3) we have with t = w µx and so Notice that the hypothesis d > 2c rules out the case We then obtain two possible values for x given by with t 1 , t 2 as in (9).If y ̸ = 0 then, from (4), γλ − γy − δw = 0, i.e.
Finally, if z ̸ = 0 then We can now write down a complete list of equilibrium points, in the half-space x ≤ w, when w > 0.
(ii.1) First, we observe that if x = y = z = 0, then from (6) we obtain w = 0, hence here we exclude this case, because we have already dealt with it (see above the case w = 0).

(ii.3)
and we obtain the equilibrium point We then obtain two equilibrium points (ii.5)If x > 0 and y = z = 0, from (6) we obtain w = 0 and we exclude this case because we are assuming x ≤ w.
Also, x is as in (10), and we obtain the two equilibrium points (ii.8)If x > 0, y > 0, z > 0, then λ − δ γ w and arguing as in [ii.7] we obtain 4 equilibrium points with w j as in (12) and t i as in (9).
In the two next sections we are going to discuss the stability of these equilibrium points.In the case w = 0, we will obtain a result independent of δ.For all the other equilibria above, we will study the stability as δ → 0 + .

Instability in the Case w = 0
In this section we study the stability of the two families of points ( 7) and ( 8).We will obtain that they all are unstable.By this we mean that for each of these points there is a neighborhood such that the trajectories starting in this neighborhood exit from it.Proposition 7. Let Z ≥ 0 and P = (0, 0, Z, 0).Then there is ϵ > 0 such that the trajectories starting in B ϵ (P) ∩ C 0 exit from B ϵ (P).
Proof.Let ϵ > 0 be such that ϵ < λδ 2(δ+σ) .This implies Let B ϵ (P) be the ball centered at P with radius ϵ and let Q = (x q , y q , z q , w q ) ∈ B ϵ (P) ∩ C 0 .As in the previous proposition, we argue by contradiction.Let u(t, Q) be the trajectory starting from Q.If it were u(t, Q) ∈ B ϵ (P) for any t ≥ 0 we would have y(t) > λ − ϵ and w(t) < ϵ for every t ≥ 0, hence Therefore, we obtain w(t) ≥ w 0 e δλ 2 t , for every t ≥ 0 and thus w → +∞, contradicting the hypothesis u(t, Q) ∈ B ϵ (P) for every t ≥ 0. The contradiction proves that there exists t > 0 such that u(t, Q) / ∈ B ϵ (P) and since this holds for every Q ∈ B ϵ (P) ∩ C 0 , our thesis is proven.

Stability of Equilibria in the Case w ̸ = 0
In this section we study the stability of equilibrium points with x ≤ w and w ̸ = 0.These are the equilibria listed above from (ii.2)-(ii.8),recalling that we have ruled out the cases (ii.1) and (ii.5).We will limit our study to what happens as δ → 0 + , assuming that the other parameters remain fixed.We will analyze the Jacobian matrix at equilibria, and this is not possible if an equilibrium lies on the hyperplane x = w, because there the vector field F is not differentiable.So let us prove, as first thing, the following proposition.Proposition 9.If P = (x 0 , y 0 , z 0 , w 0 ) ∈ C 0 is an equilibrium point and w 0 > 0, then x 0 ̸ = w 0 .
Proof.We argue by contradiction.If P is an equilibrium point and x 0 = w 0 > 0, then by (7) we obtain ≤ 1.Our hypotheses on µ rule out both these possibilities.
To study the stability of the problem, we introduce the Jacobian matrix of the differential model (1).Recall that we are studying the case x ≤ w, w > 0, and by the previous proposition we can assume x < w.We obtain As we can see, 2 is always an eigenvalue of the Jacobian calculated at the critical points.In cases from (ii.2)-(ii.4), as x = 0, we obtain ρ 1 = −c < 0. In cases (ii.6)-(ii.8)we note that where t = w/µx, so since we are in the case in which d t and when we have Thus, in cases from (ii.6)-(ii.8), the points corresponding to t 2 are always unstable, for any value of the parameters, while for equilibrium points corresponding to t = t 1 a more detailed analysis will be carried on.
Let us now study the various cases previously determined.
(ii.2) As we just said, an eigenvalue of the Jacobian matrix J(P i ) is given by ρ 1 = −c < 0.
So, we have to study the eigenvalues of the submatrix We note that we have the eigenvalue ρ 2 = γλ − δw and, for w 1 , w 2 as in (12), with k ≥ 2δ, δ → 0 + , we obtain the following expansions We need to distinguish the two cases P 1 = 0, 0, σ δ w 1 , w 1 and Because of the continuous dependence of the eigenvalues on the coefficients of the matrix, we can conclude that the point P 1 is unstable as δ → 0 + , because J(P 1 ) has the eigenvalue ρ 2 = γλ − δw 1 > 0, as δw 1 → 0 + .We have then proved the following proposition.
As to P 2 we notice that δw 2 → k for δ → 0 + , and therefore we must take in account the sign of γλ − k.If γλ − k > 0 we have that the matrix J(P 2 ) has an eigenvalue ρ 2 = ρ 2 (δ), which, by continuous dependence of the eigenvalues of a matrix on the entries, tends to the value γλ − k > 0 as δ → 0 + .We then obtain, in this case, the instability for the point P 2 as δ → 0 + .Let us now assume γλ − k < 0. We consider the submatrix Ĵ(P 2 ), defined as follows We easily obtain that, as δ → 0 + , the trace of the matrix Ĵ(P 2 ) becomes negative and its determinant positive.This implies that the eigenvalues of Ĵ(P 2 ) have strictly negative real parts.Now, the eigenvalues of the matrix J(P 2 ) are given by ρ 1 = −c < 0, ρ 2 = γλ − k and by the two eigenvalues of Ĵ(P 2 ).Hence, we have proved the following proposition.
(ii.3)We refer to P as (15).We consider J(P) for δ → 0 + and we obtain With some lengthy but standard computations we obtain the following eigenvalues As a consequence, we have the following result.
(ii.4)The Jacobian matrix calculated at the points P i , i = 1, 2 is as follows We use the Routh-Hurwitz criterion (see [11]) .As ρ 1 = −c is an eigenvalue, we can write the characteristic polynomial of J(P i ), i = 1, 2, in the following form According to the Routh-Hurwitz criterion.the necessary and sufficient conditions for all roots of p 3 (ρ) to have a strictly negative real part are the following Let us study the stability of the point P 1 obtained by setting w = w 1 as in (12).After some lengthy but standard computations, we obtain Hence, we obtain instability when kλ − σ < 0, stability when kλ − σ > 0. However, the asymptotic development of the z component of P 1 gives hence, in the case of stability, the equilibrium P 1 is not in C, as δ → 0. We conclude that the following proposition holds true.
We now study the point P 2 with w = w 2 as in (12).After some computations, we obtain the following results for the coefficients of p 3 : Hence, by applying the Routh-Hurwitz criterion again, we obtain that, as δ → 0 + , the critical point P 2 is stable when γλ − k > 0, unstable when γλ − k < 0. But in P 2 we have Hence, in the instability case γλ − k < 0, we have y < 0 as δ → 0 + , and P 2 is not in C. We have then proved the following proposition.
Proposition 14.The point w 2 as in (12), is an asymptotically stable equilibrium point if γλ − k > 0 as δ → 0 + .When γλ − k < 0 P 2 is an unstable equilibrium point of the ODE system (1), but it is not in C, as δ → 0 + .
(ii.6)In this case we have 4 critical points P i,j = P(t i , w j ), i, j = 1, 2 as in (17) with t i , w j given in ( 9) and ( 12), respectively.For sake of simplicity, we report the expression of the Jacobian matrix calculated at these points.
Based on what is established in (20), the following proposition holds.
Let us now begin the study in the case where ρ 1 (t 1 ) < 0 as seen in ( 21).Let us then analyze the eigenvalue ρ 2 = γλ − δw.We have Therefore, the Jacobian calculated at point P 1,1 has a positive eigenvalue for δ → 0 + and thus the equilibrium P 1,1 is unstable.As for point P 1,2 , we must distinguish between two cases.If γλ − k > 0 then this equilibrium will also be unstable, so let us see what happens if γλ − k < 0. In this case the eigenvalues of J(P 1,2 ) are as follows.A first one is given by ρ 1 (t 1 ) as in (21), a second one is ρ 2 = γλ − k + O(δ 2 ) < 0. The last two are the eigenvalues of the following matrix Ĵ(P 1,2 ): As δ → 0 + , this matrix has strictly negative trace and strictly positive determinant, so its eigenvalues have strictly negative real parts.Thus we have proved that, in this case, the point P 1,2 is a stable equilibrium.Let us summarize everything in the following proposition.

Let us then analyze the behavior of the Jacobian matrix calculated at point
, γλσ γσ + δ 2 , 0, δγλ γσ + δ 2 .In this case, the Jacobian matrix is almost the same one computed at point P in case (ii.3).Indeed, the two points differ only for their first entry, and the two Jacobian matrices only for the entries j 1,1 (first row, first column).Hence, the two matrices have different in the present case), but both are negative, while the other three eigenvalues are the same as in the case (ii.3).Using the results we have obtained in that case, we obtain the following proposition Proposition 18.The point P 1 = δγλ µt 1 (γσ + δ 2 ) , γλσ γσ + δ 2 , 0, γλδ γσ + δ 2 is an unstable equilibrium if kλ − σ > 0 and is an asymptotically stable equilibrium if kλ − σ < 0 as δ → 0 + .

(ii.8)
In this case we have 4 critical points P i,j = P(t i , w j ), i, j = 1, 2 as in (19) with t i , w j given in ( 9) and ( 12), respectively.Based on what is established in (20), the following proposition holds.

Proposition 19. The points P
k ≥ 2δ, d > 2c with t 2 as in (9) and w j as in (12), are unstable equilibria.
Regarding the points P 1,j (j = 1, 2) we argue as we have just done for the previous case (ii.7).Indeed, in this case, the Jacobian matrix is almost the same one computed at point P 2 in case (ii.4): the two points differ only for their first entry, and the two Jacobian matrices only for the entries j 1,1 , hence the two matrices have different ρ 1 eigenvalues, which are both negative, while the other three eigenvalues are the same in the two cases.Using the arguments and the results of case (ii.4), we obtain the following proposition Proposition 20.The point P 1,1 = w 1 µt 1 , λ − δ γ w 1 , −λ + δ γ + σ δ w 1 , w 1 , with w 1 as in (12), is an unstable equilibrium as δ → 0 + if kλ − σ < 0. If kλ − σ > 0 the point P 1,1 is an asymptotically stable equilibrium of the ODE system (1), but it is not in C, 12) is an asymptotically stable equilibrium if γλ − k > 0 as δ → 0 + .When γλ − k < 0, P 1,2 is an unstable equilibrium of ODE system (1), but it is not in C, as δ → 0 + .

Fixed Points and Stability If x > w
We are now going to look for the equilibria, and their stability, in the case x > w.To compute the equilibria, we need to consider the system given by Equation (3) and by the following equations y(γλ − γy − δx) = 0, (28) x[δy We note that x cannot be zero since x > w ≥ 0. Therefore, we will always have the relation with t i given by ( 9).Moreover, it cannot be y = z = 0 because this would lead, from (30), to w = 0 and, as just stated, x = 0. We distinguish the following cases: (iii.1) y = 0 and z ̸ = 0. From (30) we obtain z = σ δ w.From (29), we obtain Let us set w i,j = w j (t i ), with i, j = 1, 2, and t i as in (9).We have the expression for x from (31).Hence, we obtain the following equilibria (iii.2) z = 0 and y ̸ = 0. From ( 30) and (31), we have y = σ δ w = σt i µx/δ and, substituting in (28), this gives us the expression for x Then we have two equilibrium points (iii.3) y ̸ = 0, z ̸ = 0. From (29) we obtain the expression (32) for w and from (31) the expression for x.Then, (28) gives us the expression for y and finally (30) the value of z.So we have four equilibrium points Next, we will study the stability of these equilibria.For this, we compute the Jacobian matrix of the differential model ( 1) with x > w.Recalling that c = d t t 2 + 1 and considering (31) we obtain µt .Since t 1 ∈ (0, 1) and t 2 > 1 we have that if t = t 2 then j 1,1 > 0, j 1,4 < 0 and vice versa for t = t 1 .
(iii.1)As y = 0, in this case the entry j 2,2 = γλ − δx is an eigenvalue, with x = x i,j = w i,j /µt i .Let us call ρ 1 this eigenvalue.Considering the asymptotic for δ → 0 + we obtain and we can state Therefore, the Jacobian matrix calculated at the points P i,1 for δ → 0 + has at least one positive eigenvalue, and therefore the equilibria P i,1 are unstable.We thus state the following proposition.
The behavior of the points P i,2 depends instead on the sign of γλ − k.If γλ − k > 0 these equilibria are unstable.We then study the behavior of such points under the assumption γλ − k < 0. We examine the asymptotic behavior for δ → 0 + .The asymptotic of the Jacobian matrix J(P i,2 ) becomes We use the Routh-Hurwitz criterion.We introduce the characteristic polynomial of J(P i,2 ) Through some lengthy but standard computations we obtain 1) Let us recall that the necessary and sufficient condition for the roots of the polynomial p 3 (ρ) = ρ 3 + a 1 ρ 2 + a 2 ρ + a 3 to have strictly negative real parts is that a 1 , a 3 > 0 e a 1 a 2 − a 3 > 0. We recall that t 1 ∈ (0, 1) and t 2 > 1 as given in (9).Therefore, a 2 < 0 and a 3 < 0 if t = t 2 making P 2,2 unstable.Conversely, a 2 > 0 and We can therefore summarize in the following proposition.µt 1 , 0, σ δ w 1,2 , w 1,2 , with w 1,2 given in (32) and t 1 given in (9), is unstable for δ → 0 (iii.2) In this case, we have two critical points P i , i = 1, 2 as in (35).Computing the asymptotic expansions for δ → 0 + we obtain and the Jacobian matrix calculated at these points becomes , for t = t 1 and t = t 2 .When δ = 0 + , this matrix reduces to a diagonal matrix with an eigenvalue given by We have repeatedly noted that ρ 1 > 0 if t = t 2 and ρ 1 < 0 if t = t 1 .Due to the continuity of eigenvalues with respect to the elements of the matrix, in the case t = t 2 for δ → 0 + , there is an eigenvalue of the Jacobian matrix that tends to a positive value, and therefore the critical point P 2 is unstable.If t = t 1 , we obtain instead ρ 1 < 0, hence we must apply the Routh-Hurwitz criterion to the characteristic polynomial p 4 (ρ).The criterion states, in this case, that the necessary and sufficient conditions for all roots of characteristic polynomial p 4 (ρ) = ρ 4 + a 1 ρ 3 + a 2 ρ 2 + a 3 ρ + a 4 to have a strictly negative real part are the following: After lengthy computations we obtain the following asymptotic developments of the coefficients of the polynomial p 4 (ρ) Recalling that 1 − t 2 1 > 0, we have that a 1 , a 2 , a 3 > 0 for t = t 1 .Also, we have a 4 > 0 if kλ − σ < 0. Analyzing the orders of magnitude of the coefficients, the last two conditions are definitely met when the leading terms of the coefficients of the characteristic polynomial are positive.On the other hand, p 4 (ρ) obviously has a positive root if a 4 < 0. We can therefore conclude with the following proposition , for t 2 defined in (9), is an unstable equilibrium point.The point , for t 1 defined in (9), is an asymptotically stable equilibrium point if kλ − σ < 0 and is unstable if kλ − σ > 0.
(iii.3)In this case, we have four equilibrium points varying with the values of w i,j see (32) and t i , see (9).We start with w i,1 for i = 1, 2. Considering the asymptotic expansions of the variables for δ → 0 + , we obtain Notice that, when λk − σ > 0, we obtain z i,1 < 0 as δ → 0 + , hence in this case the points P i,1 are not in C.
The developments above, when substituted into the Jacobian matrix, lead to .
We now argue as in the previous case, noticing that, when δ = 0, this matrix reduces to a diagonal matrix with an eigenvalue As we know, this value satisfies ρ 1 > 0 if t = t 2 > 1 and ρ 1 < 0 if t = t 1 ∈ (0, 1).Therefore, when t = t 2 , due to the continuity of eigenvalues with respect to the entries of the matrix, there exists at least one positive eigenvalue as δ → 0 + , and thus the equilibrium point P 2,1 is unstable.We now proceed by studying the stability of the point P 1,1 .We apply again the Routh-Hurwitz criterion, and we study the characteristic polynomial p 4 (ρ) of the matrix and its coefficients a i , i = 1, .., 4. We want to study what happens for δ → 0 + , so we look at the asymptotic behavior of the a i s.We start noticing that We have that a 1 , a 2 , a 3 > 0 for t = t 1 and a 4 > 0 for t = t 1 and kλ − σ > 0. Under these conditions, the a 1 a 2 − a 3 and a 1 a 2 a 3 − a 2 1 a 4 − a 2 3 are also positive, as δ → 0 + .We recall that the equilibria P i,1 are not in C when kλ − σ > 0. We collect our results in the following proposition: Proposition 24.P 2,1 is an unstable equilibrium.If kλ − σ > 0 the equilibrium P 1,1 is asymptotically stable but is not in C, together with P 2,1 , as δ → 0 + .If kλ − σ < 0 the two equilibria are in C and are unstable, as δ → 0 + .
We now proceed with the study of the stability of the points P i,2 .The asymptotic developments in this case are which lead to the Jacobian matrix calculated at P i,2 .
We now apply the Routh-Hurwitz criterion to the characteristic polynomial p 4 (ρ) = 1) If we assume γλ − k < 0 and t = t 1 ∈ (0, 1), we have a 4 < 0, while if t = t 2 > 1 it holds a 2 , a 3 < 0, and in both cases we obtain that there is an eigenvalue with strictly positive real part.Hence, in this case, the two equilibrium points we are dealing with are unstable.Notice also that from γλ − k < 0, we derive y i,2 < 0 as δ → 0 + , so in this case the two points are not in the positive cone.
Let us now assume γλ − k > 0. If t = t 2 > 1 we obtain a 4 < 0 and the equilibrium is unstable.If t = t 1 ∈ (0, 1), all the coefficients a i are positive.For the stability we have then to verify that, at least for δ → 0 + , it holds We have that a 1 a 2 is led by a term of the type c/δ 2 , with c > 0, while a 3 = O(δ −1 ), so the first inequality above is verified for small δs.Regarding the last condition, that is a 1 a 2 a 3 − a 2 1 a 4 − a 2 3 > 0, we notice that a 2 3 = O(δ −2 ) while a 1 a 2 a 3 = O(δ −3 ) and a 2  1 a 4 = O(δ −3 ).Therefore, the last condition is simplified by requiring a 2 a 3 − a 1 a 4 > 0.
With simple calculations, we obtain so the second condition is verified for t = t 1 and δ → 0 + .We summarize in the following proposition Proposition 25.If γλ − k < 0 the two points P i,2 are both unstable equilibria of the system (1), but are not in C, as δ → 0 + .If γλ − k > 0 and δ → 0 + , the point P 2,2 is an unstable equilibrium, while the point P 1,2 is an asymptotically stable equilibrium (and both are in C).

Simulation Results
We present the results from simulations conducted with Matlab's ode45 solver, an explicit adaptive Runge-Kutta method known for its broad applicability in solving ordinary differential equations (ODEs).While our study utilizes ode45 for its efficiency and general applicability, it is worth acknowledging the existence of specialized methods for particular types of nonlinear ODEs, see for example [12], with an iterative finite difference method for nonlinear ODEs.
In all depicted graphs, the evolution of variables is illustrated using distinct colors: the population in magenta, renewable resources in green, non-renewable resources in black, and accumulated wealth in blue.Dashed lines indicate the coordinates of the equilibrium point under consideration, with matching colors for corresponding components (i.e., magenta for the population, and so forth).
In these simulations, we aimed to test our theoretical analyses by depicting in each chart a dashed line that represents the equilibrium point and a curve that results from uniformly perturbing the initial value relative to the equilibrium point.Alongside these verifications, we analyzed the behavior of the solution curves as the depletion coefficient δ of both renewable and non-renewable resources increases, sometimes obtaining results different from those achieved for a small δ (in the literature small δs are usually considered).We also analyzed the behavior of the solutions obtained by considering different depletion factors between renewable and non-renewable resources.Finally, we tested the stability of some equilibrium points supporting positive values for all the variables.On the x-axis, we use a generic unit of time that does not correspond to specific intervals such as days, months, or years, and does not necessarily span long periods.Indeed, we suspend our simulation once the expected behavior becomes clear.On the y-axis, the label "resources" refers to both renewable and non-renewable resources.This simplification was dictated by space constraints.
In these experiments, we kept certain parameter values constant, setting them to the values generally found in the literature (see [6,9]).However, we deliberately altered other parameters to study their behavior.In particular, if not stated otherwise in the comment of the individual figure, we consider the nature carrying capacity λ = 100, the nature regeneration factor γ = 0.01, the wealth consumption σ = 5, the factor that regulates the possible increase of non-renewable resources k = 1, and, finally, the factors that govern the dynamics of populations c = 1, d = 3, µ = 2. Our simulations, on the other hand, were conducted by varying the parameter δ.
In analyzing the correctness of our theoretical studies, we started from the case x ≤ w, that is, when wealth is greater than the current population.When w = 0 the critical points are given by two families of points that we have denoted by Q and R, representing the desert state and the nature state, respectively.These are two families of unstable points that in Figure 1 we have represented in the case of δ = 3, Z = 1.To achieve better chart readability with variable values of the same order of magnitude, we set λ = 10, γ = 0.5.We then analyzed the situation near the critical points of case (ii.2) both to verify the correctness of our analysis and to observe the potential scenarios around such points.The theoretical results obtained for small δ predict instability at point P 1 , whereas for point P 2 , we observe instability when γλ − k > 0 and stability otherwise.We have depicted in Figure 2a the scenario that emerges from a situation close to point P 1 = (0, 0, 9.4054, 0.9387) with δ = 0.49 < k/2.The eigenvalues of J(P 1 ) are ρ k = 0.0294, −4.7227, −1.0, 0.5316, indicating an unstable situation according to theory.This behavior was confirmed by various numerical simulations for both a small δ and a potentially large δ value.In Figure 2b, we present the expected situation following a perturbation of point P 2 when γλ − k = −0.9< 0 with λ = 10, that is, in a situation of stability.Thus, we identify the critical point P 2 = (0, 0, 10.6748, 1.0653) with eigenvalues of J(P 2 ) as ρ k = − 0.0338, −5.2929, −1.0000, −0.4316, stable according to the theory.We also examined the behavior of P 2 with a large δ value.For example, in the case δ = 5, by setting k = 2δ + 0.5 = 10.5, and thus γλ − k = −9.5, we obtained P 2 = (0, 0, 1.3702, 1.3702) and the following eigenvalues for J(P 2 ): ρ k = −3.4254± 1.6053i, −1.0000, −5.8508.The numerical results align with our predictions both when γλ − k < 0 and γλ − k > 0. We have now included an example of our simulation for large values of δ.This is just one among many tests we conducted without complete theoretical support.Because of this, we are eager to explore this area further in future research.In Figure 3, we present the analysis of point P in case (ii.3) for large values of δ.We recall that in case (ii.3), there exists a single stable critical point when kλ − σ < 0 and it becomes unstable if the condition is reversed.This analysis has been validated through numerous simulations with small δ.It should be noted that, in the case of stability, the basin of attraction of the point is very small and the perturbation applied to P had to be minimized to ensure convergence to this equilibrium point.With the standard perturbation ϵ = 10 −1 , the system tends to converge towards point P 1 of case (ii.7), which only differs from P in the population values.In Figure 3, we detail the analysis of point P = (0, 0.0500, 0, 0.1000) for δ = 10.In Figure 3a, the outcomes for kλ − σ > 0, specifically kλ − σ = 95, are shown.Here, P is established as a stable point, irrespective of kλ − σ values.However, it is crucial to highlight that the attraction basin of this point is exceedingly narrow.Similar to the previously mentioned scenario for small δ, it was necessary to reduce the standard perturbation to ϵ = 10 −3 to ensure convergence to point P instead of point P 1 = (0.1308, 0.00500, 0, 0.01000) of case (ii.7), which remains stable under identical convergence conditions for δ → 0 + .Given that P's components are near zero, a perturbation on the order of 10 −3 is deemed reasonable.The eigenvalues of matrix J(P) are ρ k = −0.2501± 0.6612i, −1.0000, −0.9896.Figure 3b explores the scenario with kλ − σ = −2, achieved with k = 0.3, λ = 10.The eigenvalues of J(P) are −0.0250± 0.0661i, −1.0000, −0.0999.Stability is observed over very long intervals.Similar to the previous cases, the stable point exhibits a notably small basin of attraction.In Figure 3b, we illustrate the scenario achieved with a perturbation ϵ = 10 −2 that leads us to convergence towards the point P 1 = (0.0131, 0.00500, 0, 0.01000) of case (ii.7).We now report some experiments conducted in the case x > w.Specifically, in Figure 4, we analyze, as δ varies, the behavior of solutions starting near the equilibrium point indicated as P 12 in the scenario (iii.1),predicted to be unstable if γλ − k > 0. In our case, we set γ = 0.1 so that γλ − k = 9 > 0.
In the case of γλ − k < 0, for which we have stability for small values of δ, we observed a situation of instability by pushing δ to the maximum value allowed while keeping k = 1 fixed.In Figure 5, we show two examples from our work on how solutions behave when the depletion rates, δ 1 and δ 2 for renewable and non-renewable resources, are different.This part of our study adds an interesting twist, showing that these different rates can change how the system behaves.We think this is a good area to explore more in future research.Specifically, in Figure 5a, we analyzed a situation similar to the one observed in Figure 3a, for the case (ii.3) with kλ − σ = 95 being unstable according to theoretical results for δ → 0 + , but exhibiting behavior that contradicts this for large δ.The point P under consideration depends only on δ 1 , but it would seem that the behavior of the solutions is dictated by δ 2 .Indeed, in Figure 5a, we report the situation obtained by setting δ 1 = 10, resulting in point P = (0, 0.0500, 0, 0.1000), and δ 2 = 0.01 .The simulation provides a scenario of unstability with eigenvalues given by ρ k = −0.2501± 0.6612i, −1.0000, 0.0089.We have also analyzed the situation with δ 1 = 0.01, δ 2 = 10 in a scenario of stability, with the eigenvalues of the Jacobian matrix in P = (0, 99.8004, 0, 0.1996) given by ρ k = − 0.9980 ± 0.0446i, −1.0000, −1.9577.Therefore, the point is stable, in accordance with the behavior observed for large values of δ.The chart is not significant as it shows a rapid convergence to the indicated point and is therefore not included here.
In Figure 5b, we instead report the scenario obtained for point P 1 in case (iii.2),again with δ 1 ̸ = δ 2 when kλ − σ = 95 > 0. It is recalled that the theoretical analysis conducted for δ 1 = δ 2 = δ predicts a situation of instability for small δ.However, numerical experiments have revealed a behavior of stability for large δ.Returning to δ 1 ̸ = δ 2 , even in this case, P 1 depends only on the value of δ 1 .In the Figure 5b we have reported the result obtained with δ 1 = 3 leading to P 1 = (0.3319, 0.4226, 0, 0.2536) and δ 2 = 0.1.This is a situation of instability, and the eigenvalues of the Jacobian matrix are ρ k = −0.1567± 0.7521i, −2.0957, 0.0444.The presence of a small δ 2 thus leads us to achieve what is predicted by the theory for δ → 0 + .We do not report the graph of the result obtained by inverting the values of δ 1 and δ 2 , that is, by setting δ 1 = 0.1 leading to P 1 = (2.0748,79.2516, 0, 1.5850) and δ 2 = 3.In this case, we achieved a situation of stability with the eigenvalues of the Jacobian matrix given by ρ k = −0.7603± 0.4076i, −10.3915, −5.4577.The graph turned out to be of little significance as the stability situation was reached quickly.We therefore conclude by saying that even in the case (iii.2), as in the case of (ii.3), it would seem that δ 2 controls the behavior of the equilibrium point whose coordinates depend only on δ 1 .
We would have liked to further analyze the scenario presented in (ii.4), because the critical point depends on both δ 1 and δ 2 , but studies for significantly different values of the two parameters are not possible because, as they increase, the point exits the cone.
A state of stability with a positive population is depicted in Figure 6.We have analyzed the scenario described in cases (ii.8) and (iii.3)under the predicted conditions of stability as δ → 0 + .In both instances, we set δ = 0.3.Figure 6a