The world faces an alarming issue today due to the depletion of forest resources caused by deforestation, fires, illegal logging, and other factors. Many countries will lose their remaining forests by 2030 if this trend continues, according to a recent report [1]. Urgent action is needed to address this challenge, including better coordination and control of the timber industry and communities that depend on the forests [2]. Mathematics provides some powerful tools to tackle such problems with the help of differential equations which can offer a way to solve dynamical systems making them essential to science, engineering and humanity. Some studies have used the mathematical modeling of forest depletion and suggested solutions using various numerical and analytical methods. Gompil et al. [3] proposed numerical and simulated results for a forest depletion model, while Eswari et al. [4] examined the homotopy perturbation method (HPM) to solve the mathematical model for the depletion of forest resources. Nugraheni et al. [5] proposed stability analysis and numerical simulations for a mangrove forest resource dynamical model. Didiharyono and Kasse [6] studied the stability of a mathematical model for deforestation and presented numerical simulations of the system. All these studies offer useful insights into the dynamics of forest resources and propose possible solutions to handle this critical global problem. This paper concentrates on the study of the depletion of forest resources, employing a mathematical model suggested by Misra, Lata, and Shukla [7]. This mathematical model consists of the cumulative density of forest resources, the density of the population, and population pressure, represented by the variables B, N, and P, respectively. In this model, the connection between forest and population density is considered as a prey-predator logistic model. The forest density decreases as housing and development increase, impacting its growth rate. Population pressure growth is proportional to population density in the model [7]. The authors have investigated existence and uniqueness of the global positive solution and provided numerical simulations to study this model. The cumulative density of forests and population size, are modelled using comprehensive equations with dynamic relations similar to a prey-predator system. The model signifies the depletion of forest resources provoked by population growth, reduction of forest areas for expansion purposes and the depletion by the pressure of the population. In addition, the model considers that the increase in population pressure is proportional to population density. This model consists of dimensionless differential equations. The suggested model [7] can be represented as: d B d t = s B − h B 2 − α B N − λ 2 B 2 P , d N d t = r N − j N 2 + π α B N , d P d t = λ N − λ 0 P , (1) where B (0) ≥ 0, N (0) ≥ 0, P (0) ≥ 0 and we define variables and constant coefficients of this model in the following table as.

Values for the parameters and coefficients are considered, s = 0.8, s ₀ = 0.2, L = 50, α = 0.0001, λ = 0.2, λ ₀ = 0.1, r = 0.2, r ₀ = 0.1, K = 100, π = 0.004, λ ₂ = 0.0002, h = s 0 L , j = r 0 K and initial conditions n ₁ = B (0) = 30, n ₂ = N (0) = 35, n ₃ = P (0) = 1 as given in (Misra et al., 2014).

The rate of forest depletion is alarming, mainly driven by illegal logging and land clearing activities. This trend poses a serious threat to our ecosystem and immediate action is needed to mitigate its impact. Without decisive measures, the depletion of forest resources will have long-lasting consequences on the environment and our wellbeing. It is imperative to find sustainable solutions and enforce regulations to curb the rampant destruction of forests and preserve them for future generations.

Many dynamical problems in science and engineering cannot be solved analytically (exactly) and can be approximated numerically. There is another technique named as series solution (semi-analytical techniques) which is more closer to analytical results. For this purpose a wide range of methods have been developed to find approximate solutions that are as close as possible to the exact solutions. Among these methods are the Taylor series method [8], which approximates functions as power series; the Picard method [9], which iteratively computes solutions from initial conditions; the Adomian decomposition method [10], which decomposes a differential equation into simpler sub problems; the variational iteration method [11], which uses Lagrange multipliers to optimize solutions; and the homotopy perturbation method [12–14,14,15,17–19], which constructs a homotopy that gradually deforms the problem into a simpler one while adding a perturbation term to the solution. These methods have been applied to a wide range of problems in physics, engineering, various fields and have proven to be highly effective in providing accurate approximations to complex dynamical systems.

Ji Huan He, a mathematician from China proposed a novel semi-analytical method based on homotopy and perturbation techniques in 1999, which was named the homotopy perturbation method (HPM) [12]. He improved and extended the HPM to solve a wide range of problems. In 2004, He used the HPM to non-linear oscillators and asymptotic [13,14]. In 2005, the HPM was applied to solve non-linear wave equations and problems related to limit cycle and bifurcation of non-linear systems [15,16]. In 2008, He employed the HPM to solve boundary value problems [20]. In 2007, Javidi and Golbabai used a revised version of the HPM to solve non-linear Fredholm integral equations [21].Recently, HPM with small variations has been applied to study fractal duffing oscillator problems under arbitrary conditions [22], modified HPM for nonlinear oscillators Anjum and He [23], attachment oscillator arising in nanotechnology [24], conservative nonlinear oscillators [25], non-linear oscillator problems in a fractal space [26] and HPM including Aboodh transformation to solve fractional calculus Tao et al. [27], vibrating magnetic inverted pendulum Moatimid et al. [28], Symmetry-breaking and pull-down motion for the helmholtz–duffing oscillator Niu et al. [29], nonlinear fractional Drinfeld–Sokolov–Wilson Equation Nadeem and Alsayaad [30], trajectory analysis of a zero-pitch-angle e-Sail Niccolai et al. [31], natural convection between two concentric horizontal circular cylinders Abdulameer and Ali Al-Saif [32], nonlocal initial-boundary value problems for parabolic and hyperbolic Al-Hayani and Younis [33], multi-step iterative methods for solving nonlinear equations Saeed et al. [34], telegraph equation Moazzzam et al. [35], triangular linear diophantine fuzzy system of equations Shams et al. [36], condensing coagulation model and Lifshitz-Slyzov equation Arora et al. [37], singular nonlinear system of boundary value problems Pathak et al. [38], rikitake-yype system Ene and Pop [39], heat and mass transfer with 2D unsteady squeezing viscous flow problem Abdul-Ameer and Ali Al-Saif [40], variable Speed Wind Turbine Control Shalbafian and Ganjefar [41], radial thrust problem Niccolai et al. [42], special third grade fluid flow with viscous dissipation effect over a stretching sheet Swain et al. [43], and the frequency–amplitude relationship of a nonlinear oscillator with cubic and quintic nonlinearities He et al. [44]. The HPM has become a widely-used technique to solve a large variety of problems in different fields and many research papers have been published each year using this method as evidenced by a simple search on Google Scholar.

In this paper, we provide an approximate solution of model 1) by using the homotopy perturbation method. The interesting feature of HPM is that it provides the best approximate solution by taking a few numbers of perturbation terms.

Consider a non-linear differential equation D μ − g τ = 0 , τ ∈ ℧ (2) subject to the boundary condition β μ , ∂ μ ∂ τ = 0 , τ ∈ Γ (3) where D is a differential operator, β is boundary operator, Γ is the boundary of the domain ℧ and g(τ) is an unknown function. The D, generally consist on two parts, linear and non-linear part, represented as L and N respectively. Therefore, 2) can be written as follows L μ + N μ − g τ = 0 , (4) using homotopy method, by taking an embedding parameter q one can construct a homotopy v (τ, q): ℧ × [0, 1] → R for Eq. 4 which satisfies H w , q = 1 − q L w − L μ 0 + q L w + N w − g τ = 0 , (5)

and it is equivalent to H w , q = L w − L μ 0 + q L μ 0 + q N w − g τ = 0 , (6)

where q ∈ [0, 1] is an embedding parameter, μ ₀ is an initial guess approximation of Eq. 6 which satisfies the initial (or boundary) conditions. It can be written as follows. q = 0 ,           H w , 0 = L w − L μ 0 , (7) q = 1 ,           H w , 1 = L w + N w − g τ . (8)

We suppose the solution in the form of power series for Eq. 5 by taking an embedding parameter q w = w 0 + q w 1 + q 2 w 2 + q 3 w 3 + ⋯ (9)

The approximate solution of Eq. 2 can be obtained by setting q = 1, μ = lim q → 1 w = w 0 + w 1 + w 2 + w 3 + ⋯ (10)

The convergence of (Eq. 10) has been proved in [12]. The series is convergent for most cases, however, the convergent rate depends upon the nonlinear operator N(w). Furthermore He suggested the following conditions.

1. The second derivative of nonlinear operator N(w) with respect to w must be small, because the parameter q may be relatively large, i.e., q → 1.

Now we apply HPM on our model, Eq. 1 of depletion of forest resources (non-linear system of differential equations) as 1 − q u ′ − B 0 ′ + q u ′ − s u + h u 2 + α u v + λ 2 u 2 w = 0 , 1 − q v ′ − N 0 ′ + q v ′ − r v + j v 2 − π α u v = 0 , 1 − q w ′ − P 0 ′ + q w ′ − λ v + λ 0 w = 0 . (11)

The initial guesses for (11) are constant as defined in [7]. u 0 t = B 0 t = B 0 = n 1 v 0 t = N 0 t = N 0 = n 2 w 0 t = P 0 t = P 0 = n 3 (12)

and we assume the solution of (11) as, u = u 0 + q u 1 + q 2 u 2 + q 3 u 3 + ⋯ , v = v 0 + q v 1 + q 2 v 2 + q 3 v 3 + ⋯ , w = w 0 + q w 1 + q 2 w 2 + q 3 w 3 + ⋯ , (13)

by substituting Eq. 13 in Eq. 11 and collecting the terms of powers of q, we obtain q 0 : u 0 ′ = 0 ,    u 0 0 = n 1 , v 0 ′ = 0 ,    v 0 0 = n 2 , w 0 ′ = 0 ,    w 0 0 = n 3 . (14) q 1 : u 1 ′ + u 0 α v 0 − s + u 0 2 h + λ 2 w 0 = 0 ,     u 1 0 = 0 , v 1 ′ − r + α π u 0 v 0 + j v 0 2 = 0 ,    v 1 0 = 0 , w 1 ′ − λ v 0 + λ 0 w 0 = 0 ,    w 1 0 = 0 . (15) q 2 : u 2 ′ + α u 0 v 1 + u 1 α v 0 + 2 u 0 h + λ 2 w 0 − s + λ 2 u 0 2 w 1 ,    u 2 0 = 0 , v 2 ′ − α π u 1 w 0 − r + α π u 0 − 2 j v 0 v 1 = 0 ,    v 2 0 = 0 , w 2 ′ − λ v 1 + λ 0 w 1 = 0 ,    w 2 0 = 0 , (16) q 3 : u 3 ′ + α u 0 v 2 + u 1 2 h + λ 2 w 0 + u 2 α v 0 + 2 u 0 h + λ 2 w 0 + u 1 α v 1 + 2 λ 2 u 0 w 1 + 2 λ 2 u 0 2 w 2 = 0 ,      u 3 0 = 0 , v 3 ′ − α π u 2 v 0 − α π u 1 v 1 + j v 1 2 − r v 2 − α π u 0 v 2 + 2 j v 0 v 2 = 0 , v 3 0 = 0 , w 3 − λ v 2 + λ 0 w 2 = 0 ,      w 3 0 = 0 . (17) q 4 : u 4 ′ + α u 2 v 1 + α u 0 v 3 + u 3 − s + α v 0 + 2 u 0 h + λ 2 w 0 + λ 2 u 1 2 w 1 + λ 2 u 0 u 2 w 1 + u 1 α v 2 + 2 u 2 h + λ 2 w 0 + 2 λ 2 u 0 w 2 + λ 2 u 0 2 w 4 = 0 ,      u 4 0 = 0 , v 4 ′ − α π u 3 v 0 − α π u 2 v 1 − α π u 1 v 2 + 2 j v 1 v 2 − r v 3 − α π u 0 v 3 + 2 j v 0 v 3 = 0 ,      v 4 0 = 0 , w 4 − λ v 3 + λ 0 w 4 = 0 ,      w 4 0 = 0 . (18) . . .

Now considering s = 0.8, s ₀ = 0.2, L = 50, α = 0.0001, λ = 0.2, λ ₀ = 0.1, r = 0.2, r ₀ = 0.1, K = 100, π = 0.004, λ ₂ = 0.0002, h = s 0 L , j = r 0 K , n ₁ = B (0) = 30, n ₂ = N (0) = 35, n ₃ = P (0) = 1, and simplifying the equations from (Eqs 14–18) we have.

By substituting these values in Eq. 13, we have the solution of model 1) as B t = lim q → 1 u = 30 + 20.115 t + 4.84665 t 2 − 0.260167 t 3 − 0.495765 t 4 − 0.12174 t 5 + ⋯ , (19) N t = lim q → 1 v = 35 + 5.77542 t + 0.375578 t 2 + 0.00519615 t 3 − 0.000913025 t 4 − 0.00006 t 5 + ⋯ , (20) P t = lim q → 1 w = 1 + 6.9 t + 0.232542 t 2 + 0.0172871 t 3 − 0.00017237 t 4 − 0.000033 t 5 + ⋯ (21)

For α = 0.0001, α = 0.0002 and α = 0.0004, we have.

B(t) _α=0.0001 = 30 + 20.115t + 4.84665t ² − 0.260167t ³ − 0.495765t ⁴ + ⋯,

B(t) _α=0.0002 = 30 + 20.01t + 4.77438t ² − 0.274268t ³ − 0.49119t ⁴ + ⋯ and.

B(t) _α=0.0004 = 30 + 19.8t + 4.63094t ² − 0.301744t ³ − 0.481988t ⁴ + ⋯.

For λ = 0.1, λ = 0.2 and λ = 0.3, we have.

B(t) _λ=0.1 = 30 + 20.115t + 5.16165t ² + 0.0854419t ³ − 0.336328t ⁴ + ⋯,

B(t) _λ=0.2 = 30 + 20.115t + 4.84665t ² − 0.260167t ³ − 0.495765t ⁴ + ⋯ and.

B(t) _λ=0.3 = 30 + 20.115t + 4.53165t ² − 0.605776t ³ − 0.648587t ⁴ + ⋯.

For λ ₂ = 0.0001, λ ₂ = 0.0002 and λ ₂ = 0.0003, we have.

B ( t ) λ 2 = 0.0001 = 30 + 20.205 t + 5.24226 t 2 + 0.113954 t 3 − 0.334519 t 4 + ⋯ ,

B ( t ) λ 2 = 0.0002 = 30 + 20.115 t + 4.84665 t 2 − 0.260167 t 3 − 0.495765 t 4 + ⋯ and.

B ( t ) λ 2 = 0.0002 = 30 + 20.025 t + 4.45157 t 2 − 0.629906 t 3 − 0.645153 t 4 + ⋯ .

In this paper, we used the homotopy perturbation method to obtain a semi-analytical solution for the nonlinear model of the depletion of forest resources. Important characteristic of HPM is that it provides the adequate approximate series solution by taking a few number of perturbation terms which is near to analytical exact solution. Through comparison with the Runge-Kutta method, we established the effectiveness and accuracy of the proposed method. Additionally, we investigated the behaviour of the model by varying the values of the depletion rate of forest resources due to population α, the growth rate coefficient of population pressure caused by population λ, and the depletion rate of its carrying capacity due to population pressure λ ₂. The results showed that reducing these coefficients can increase the cumulative density of forest resources B(t). These findings highlight the urgent need for measures to conserve forest resources for the wellbeing of our planet. The presented model and its solution indicate the seriousness of this global issue which needs to be acted upon immediately and effectively to preserve our forest resources. This study suggests that additional investigations and research is needed to build more relevant models for assistance of forest resource experts.