Comprehensive Study of a Diabetes Mellitus Mathematical Model Using Numerical Methods with Stability and Parametric Analysis

Diabetes is sweeping the world as a silent epidemic, posing a growing threat to public health. Modeling diabetes is an effective method to monitor the increasing prevalence of diabetes and develop cost-effective strategies that control the incidence of diabetes and its complications. This paper focuses on a mathematical model known as the diabetes complication (DC) model. The DC model is analyzed using different numerical methods to monitor the diabetic population over time. This is by analyzing the model using five different numerical methods. Furthermore, the effect of the time step size and the various parameters affecting the diabetic situation is examined. The DC model is dependent on some parameters whose values play a vital role in the convergence of the model. Thus, parametric analysis was implemented and later discussed in this paper. Essentially, the Runge–Kutta (RK) method provides the highest accuracy. Moreover, Adam–Moulton’s method also provides good results. Ultimately, a comprehensive understanding of the development of diabetes complications after diagnosis is provided in this paper. The results can be used to understand how to improve the overall public health of a country, as governments ought to develop effective strategic initiatives for the screening and treatment of diabetes.


Introduction
Diabetes is sweeping the globe as a silent epidemic, posing a rising threat to public health. On a global scale, there has been a significant increase in the prevalence and incidence of diabetes. According to the most recent 2021 data from the International Diabetes Federation (IDF), 537 million adults, aged 20-79 years, are living with diabetes [1]. This number is estimated to increase to 784 million by 2045. Moreover, diabetes caused 6.7 million deaths in 2021, which is one in every five seconds [1]. Diabetes is the second most prevalent comorbidity in COVID-19, according to epidemiological research [2]. Diabetes patients are more likely to have disease severity, slower recovery, and COVID-19 consequences, including ICU admission, mechanical ventilation, and even death [3]. These issues indicate the dire need for effective intervention techniques and policies critical to the prevention of the rise in the number of individuals with diabetes [1].
Diabetes is a non-communicable, long-term illness characterized by elevated blood sugar levels. Diabetes may affect anyone, irrespective of size, age, or gender. Diabetes can be mainly categorized into two types. Type 1 diabetes is a chronic disease in which This paper focuses on a mathematical model simulating the population dynamics of diabetes and its complications. The model was introduced by Boutayeb et al. in 2004, also known as the diabetes complication (DC) model [22]. It is divided into two compartments, diabetics with complications (C) and diabetics without complications (D). In their study, two methods were used to analyze the DC diabetes model: Method I and Euler [5]. Similarly, Akinsola and Oluyo used Euler and Range-Kutta to analyze the DC model [28]. This study used the explicit Euler method, implicit Euler method, Heun's method, the 4th order Runge-Kutta method, and the Adam-Moulton method to analyze the DC model. The DC model is analyzed using different numerical methods to monitor the size of the diabetic population with and without complications over time. This will provide intriguing opportunities for developing and testing theories, estimating and regulating parameters, comprehending the dynamics of the population, and suggesting practical and effective preventive measures depending on various scenarios [19]. The effect of a growing or declining incidence of diabetes and its complications can be visualized using parameters. The model demonstrates the various strategies that may be developed by varying the parameters that characterize the incidence of diabetes, and diabetes-related comorbidities.
The contributions of this paper are summarized as follows: • Investigation of five different numerical methods to analyze the diabetes mathematical model; two have not been used before for this problem. This allows for a comprehensive comparison of the numerical methods and an understanding of the best method to analyze the biological system; • Examination of the system behavior for different values of all the parameters, which has not been addressed in the literature. Based on the results, we comment on the best values of the parameters for stability analysis; • Investigating the system for different values of diabetes incidence. The incidence of diabetes with complications is significant and of concern, necessitating early and effective strategic initiatives for screening, treatment, and management of diabetes.

Materials and Methods
Stability, analytical, and numerical methods were used to examine the model. The codes were made using MATLAB software to generate graphs and numerical data. Figure 1 illustrates the DC mathematical model proposed by Akinsola and Temitayo [29]. This is the model that was used in this paper for numerical analysis and to study the control of diabetes epidemiology over the years. diabetes diagnosis [19,26]. According to Gatwood et al., more than 30% of veterans had chronic kidney disease before being issued a diabetes diagnosis [27]. This paper focuses on a mathematical model simulating the population dynamics of diabetes and its complications. The model was introduced by Boutayeb et al. in 2004, also known as the diabetes complication (DC) model [22]. It is divided into two compartments, diabetics with complications (C) and diabetics without complications (D). In their study, two methods were used to analyze the DC diabetes model: Method I and Euler [5]. Similarly, Akinsola and Oluyo used Euler and Range-Kutta to analyze the DC model [28]. This study used the explicit Euler method, implicit Euler method, Heun's method, the 4th order Runge-Kutta method, and the Adam-Moulton method to analyze the DC model. The DC model is analyzed using different numerical methods to monitor the size of the diabetic population with and without complications over time. This will provide intriguing opportunities for developing and testing theories, estimating and regulating parameters, comprehending the dynamics of the population, and suggesting practical and effective preventive measures depending on various scenarios [19]. The effect of a growing or declining incidence of diabetes and its complications can be visualized using parameters. The model demonstrates the various strategies that may be developed by varying the parameters that characterize the incidence of diabetes, and diabetes-related comorbidities.
The contributions of this paper are summarized as follows: • Investigation of five different numerical methods to analyze the diabetes mathematical model; two have not been used before for this problem. This allows for a comprehensive comparison of the numerical methods and an understanding of the best method to analyze the biological system; • Examination of the system behavior for different values of all the parameters, which has not been addressed in the literature. Based on the results, we comment on the best values of the parameters for stability analysis; • Investigating the system for different values of diabetes incidence. The incidence of diabetes with complications is significant and of concern, necessitating early and effective strategic initiatives for screening, treatment, and management of diabetes.

Materials and Methods
Stability, analytical, and numerical methods were used to examine the model. The codes were made using MATLAB software to generate graphs and numerical data. Figure  1 illustrates the DC mathematical model proposed by Akinsola and Temitayo [29]. This is the model that was used in this paper for numerical analysis and to study the control of diabetes epidemiology over the years.  The model examines parameters such as the natural mortality rate (µ), probability of developing a complication (λ), the rate at which complications are controlled (γ), the rate at which patients with complications become severely disabled (v), and mortality rate due to complications (δ). All those parameters affect the model; however, parameters µ, γ, v, and δ are natural rates that describe death or control to patients with diabetes. The only parameter that can be changed is λ, which could be controlled by providing society with motivation to adopt good habits and develop policies such as diet, exercise, regular medical checkups, and more. Therefore, the model examines and analyzes the effect of changing this value that correlated to developing complications, to see its overall effect on the system. The effect of different values of other parameters is also explored. The model is implemented while keeping the incidence of people with diabetes constant. At the time of diagnosis, many patients already have microvascular and macrovascular complications [26], which is reflected in the model. In general, the model's independent variable is time. In addition, the dependent variables are the number of people who have diabetes without complication D(t) and the number of people who have diabetes with complication C(t). N(t) is the total number of people with diabetes.

ODEs for Modeling Diabetes
The model is explained mainly through two differential equations. One is for diabetics with complications (1) and the other is for non-complication diabetics (2).
Formulas (3) and (4) made a first-order system of differential equations where θ = γ + µ + v + δ , N (t) = dN(t) dt and N(t): is the total number of diabetics (N = C + D). The initial conditions of C(t) and N(t) are denoted as C 0 and N 0 .
Thus, the analytical solutions (5) and (6) were used to find the true value of C(t) and N(t). Hence, the true values were used for comparison and computing the percentage error.

The Parameters of the ODEs
The parameters of the ODEs play an important role in identifying major changes in diabetic patients over time, with respect to the incidence of diabetes. The critical points of N(t) and C(t) as has been found by Akinsola and Temitayo [29] to be as (15) and (16): Then the initial values, as suggested by [22], are (17) and (18): The probabilistic parameters and rates of the ODE model are defined in Table 1. These parameters were used to analyze the numerical solutions of the model. It was observed that most of the authors [19,22,28,30] define γ as the rate at which complications are cured. Similar to Akinsola and Temitayo [29], for ease of comprehension, this paper will refer γ to the rate at which complications are controlled due to the understanding that most of the complications of diabetes are chronic and not curable. The other complications are acute and mostly medical emergency cases [31].

Stability Analysis
If critical points are applied to the system, it reaches a steady state. Through that, we were able to obtain the characteristic equation of the model (21). The quadratic equation is used to find the discriminant (22) and ultimately calculate the eigenvalues (23). Furthermore, the eigenvalues can be analyzed to identify the type of critical point and type of stability of the system. If the determinant is (∆ > 0), then eigenvectors are real and negative. Furthermore, if ∆ = 0, both eigenvectors equal each other, considered real and negative. Finally, if ∆ < 0, both eigenvectors are complex conjugates with negative parts. Finally, the type of critical point and stability of the system are shown in Table 2. Using (15)- (18), the eigenvalue is equal to χ 1 = −0.156665 and χ 2 = −0.943351. Furthermore, our ∆ = 0.6189. According to Table 2, For the determinant, if (∆ > 0), then the eigenvectors are both real and negative [5]. Furthermore, Table 2 suggests that the model is asymptomatically stable. Table 2. Type of critical points and stability based on the eigenvalues and determinants [5].

Type of Critical Point Stability
Proper or improper node Unstable Proper or improper node Asymptomatically stable Asymptomatically stable χ 1 and χ 2 are the eigenvalues, λ is the probability of developing a complication, and µ is the natural mortality rate.

Numerical Methods Analysis
Numerical analysis is the process of generating numerical solutions to mathematical expressions using numerical methods [29]. It entails developing, analyzing, and implementing computer systems to solve continuous numerical problems [29]. If the numerical solution approaches the exact solution as the step size approaches zero, the numerical method is convergent [32]. Furthermore, the method must be both convergent and stable because a slight disturbance of the input data does not disrupt the convergence and causes only a minor increase in error [29,33]. The numerical analysis of the model was performed using the MATLAB environment.

Explicit Euler Method
In the explicit method, the forward finite difference O(h) was used to represent the system of ODEs by obtaining an approximate solution.
Although the explicit method is rather simple, it has the limitation of being conditionally stable. To solve this issue, the amplification factor is suggested to be equal approximately to 1 (when we compare the new value with the one before) to improve the stability of the ODEs in the explicit method and the error as well, since this method has the highest percentage of error out of the other methods. where they have a common factor of C i+1 = −αC i ∆t + C i , (similarly applied for N as well), thus, the amplification factor (G) is equal to: The discretization of Equations (3) and (4) using the Explicit Euler method results in the following numerical equations:

Implicit Euler Method
The explicit method calculates the future or the new value from the currently known system status. However, the implicit method calculates the future value from the system at present and future times. When the state of C i and N i is known, C i+1 and N i+1 can be calculated. Ultimately, the implicit method is unconditionally stable and allows the use of larger step sizes (∆t).
The discretization of Equations (3) and (4) using the implicit Euler method results in the following numerical equations:

Heun's Method
Heun's method is also another Euler method that is both explicit and implicit. Although it is also based on Euler methods, it provides higher accuracy than both. This method provides an improved accuracy compared to the last two methods due to improving the slope estimation and determination of the two derivatives, one at the beginning of the interval and one at the end. Then, these two derivatives are averaged together to obtain a better slope estimation. Ultimately, Heun's method is a modified Euler method using the predictor and corrector equations for better results. Heun's method is a 2 nd -order error O h 2 method, which yields better error percentages.
The discretization of Equations (3) and (4) using Heun's method results in the following numerical equations:

Runge-Kutta 4th Order Method
The Runge-Kutta (RK) method consists of an ODE that defines dC dt and dN dt at initial time t(0). The RK method is useful since it is able to find an unknown value of function t (time) at any C or N. The formulas (34)-(37) are used for calculating the next value for C i+1 and N i+1 from the previous value C i and N i . ∆C 1 represents the slope at the beginning, ∆C 2 represents the slope occurring at the midpoint of the interval between t and ∆C 1 , ∆C 3 represents the slope of the midpoint interval between t and ∆C 2 , ∆C 4 represents reaching the end of the interval, and similarly for ∆N 1,2,3,4 . The RK method is said to be stable if the Eigen values are real without a complex conjugate, less than zero, which is the case in this study.

Adam-Moulton Method
The Adam-Moulton method of the 4th order allows us to explicitly compute the approximate solution at an instant time from the solution in previous instants. As the initial value at the initial time step is known, then the RK method can be used to get the midpoints at the intervals (C 1 , C 2 , C 3 , and C 4 ), and apply them to the predictor and corrector formulas of Adam-Moulton. The predictor formula is used to calculate a rough approximation of C i and N i from the RK method, and then approximate the solution by using the corrector formula.
The RK method is used to find C k and N k when k = 0,1,2,3: Then the predictor and corrector equations for k = 4:

Results
We implemented the five methods and compared the results with the true values. By calculating the error, the accuracy of the proposed model for the methods used can be evaluated. As real changes in the prevalence will require years to be noted, it is recommended to begin calculating the prevalence after five years of the current time. In addition, the time step (∆t) has been chosen to be one year. Therefore, once the prevalence is noted, the government can plan properly and act quickly on implementing different measures for the sake of the population. Table 3 shows the results for C i and N i using the explicit Euler, implicit Euler, Heun's, Adam-Moulton, and RK methods from five years and above with ∆t = 1 year. Noting that the years that we chose to include are only the years 5, 10, 15, and 20, as significant changes can be noticed only every couple of years. The table shows that the RK method is more efficient than the others for this mathematical model.

Stability Analysis
The following sections will provide a comprehensive analysis of the effects of varying the time step size and the parameters of the model on the stability of the system. The changes in ∆t affect the accuracy of the results. Therefore, the lower the value of ∆t, the better the results. Table 4 compares the results of the error rate between ∆t = 1 and ∆t = 0.5, for all the methods for C i only, while Table 5 compares the results of N i . It can be seen that the error rates are much lower for ∆t = 0.5, resulting in better accuracy rates of the model. C i is the number of diabetic patients with complications, E(C i ) is the difference between C i and its true value, t is the time in years, and ∆t is the time step size. C i is the number of diabetic patients with complications, E(C i ) is the difference between C i and its true value, t is the time in years, and ∆t is the time step size.
As the stability is hugely affected by the step size (∆t), this paper aims to study the stability of all the methods used while changing the ∆t. The values of ∆t used in this paper range from 0. 1 to 1.5. Figure 2 shows the results for both C i and N i for (a) explicit Euler, (b) implicit Euler, (c) Heun's, (d) Adam-Moulton, and (e) RK methods. It can be seen that some methods are diverging, beginning from ∆t = 1.5 and above. In addition, the results with smaller ∆t are much more accurate and stable.  is the number of diabetic patients with complications, E( ) is the difference between and its true value, t is the time in years, and ∆ is the time step size. is the number of diabetic patients with complications, E( ) is the difference between and its true value, t is the time in years, and ∆ is the time step size.
As the stability is hugely affected by the step size (∆ ), this paper aims to study the stability of all the methods used while changing the ∆ . The values of ∆ used in this paper range from 0. 1 to 1.5. Figure 2 shows the results for both Ci and Ni for (a) explicit Euler, (b) implicit Euler, (c) Heun's, (d) Adam-Moulton, and (e) RK methods. It can be seen that some methods are diverging, beginning from ∆ = 1.5 and above. In addition, the results with smaller ∆ are much more accurate and stable.

The Effect of Varying the Parameters Values
The parameters of the system defined in Table 1 are specific values that ensure the stability of the system. In addition, the stability of the model using the RK method was evaluated while altering the parameters mentioned in Table 1, which are γ, λ, δ, µ, and v.

The Rate at Which Complications Are Controlled (γ)
The investigated γ values in this paper are 0.04, 0.08, 0.5, and 1.0. These values were obtained from [22]. Figure 3 shows the results of C i and N i for the RK method in terms of γ and ∆t, where (a) γ = 0.04, (b) γ = 0.08, (c) γ = 0.5, and (d) γ = 1.0. It can be seen that as larger the γ value, the larger the error rates. After γ = 1.0, the results started to diverge. (c) Heun's method; (d) RK method; (e) Adam's method, where is the number of diabetic patients with complications and is its true value, is the total number of people with diabetes and is its true value, and d is the time step size.

The Effect of Varying the Parameters Values
The parameters of the system defined in Table 1 are specific values that ensure the stability of the system. In addition, the stability of the model using the RK method was evaluated while altering the parameters mentioned in Table 1, which are , λ, , , and .

The Rate at Which Complications Are Controlled (γ)
The investigated γ values in this paper are 0.04, 0.08, 0.5, and 1.0. These values were obtained from [22]. Figure 3 shows the results of Ci and Ni for the RK method in terms of γ and ∆ , where (a) γ = 0.04, (b) γ = 0.08, (c) γ = 0.5, and (d) γ = 1.0. It can be seen that as larger the γ value, the larger the error rates. After γ = 1.0, the results started to diverge.  the total number of people with diabetes and is its true value, and d is the time step size.

Probability of Developing a Complication (λ)
The investigated λ values in this paper are 0.04, 0.66, 0.85, and 1.2. Figure 4 shows the results of Ci and Ni for the RK method in terms of λ and ∆ , where (a) λ = 0.04, (b) λ = 0.66, and (c) λ = 0.85.   is the number of diabetic patients with complications and is its true value, is the total num ber of people with diabetes and is its true value, and d is the time step size.
The Mortality Rate Due to Complications ( ) The investigated values in this paper are 0.02, 0.05, 0.1, and 0.35. Figure

Discussion
The mathematical modeling of diabetes aids in the analysis and interpretation of population changes. In previous studies, some papers have solved the diabetic model using parameters related to the specific country under study. This paper adds to the literature by performing a comprehensive evaluation of the effect of the time step and all the possible values of relevant parameters using multiple methods that can be applied to any country.
In this study, two ODEs were to analyze the development of diabetic patient numbers with complications. The ODEs were solved using five numerical methods, namely, explicit Euler, implicit Euler, Heun's, RK (4th order), and Adam-Moulton (4th order). The implicit, explicit, and RK methods were used in previous studies, and in this paper, Heun's, and Adam-Moulton's methods are presented for the first time to study and analyze the prevalence of diabetes. After discretization and implementation of all methods, the RK method resulted in the best accuracy rates, with a low 1 × 10 −6 % error rate. RK method is compared to other implemented methods in this paper due to its superior performance. particularly, the performances of the RK and the Adam-Moulton methods are always compared for first-order ODEs. According to Gofe and Gebregiorgis [30], the Adam-Moulton method is considered the most efficient due to the predictor and modifier equations in the model. However, in terms of accuracy, there is no generalization on the best method, and it depends solely on the ODEs used. As a result, it is recommended to observe the performance of each method based on the performance of ODE after doing the stability analysis. Moreover, in the literature, the implicit Euler method showed a smaller percentage error than the explicit Euler. This paper highlights similar results [23]. Also, the Implicit method proved more accurate than Heun's method for the problem at hand.
It can be noticed that the model is stable using implicit Euler, Heun's, and RK's methods, as they perform and calculate implicitly. The explicit Euler and Adam-Moulton's methods are stable for ∆t less than 1.0 and begin to become unstable for ∆t values larger than one year, as shown in Figure 2; this is because both work explicitly. When ∆t = 1.5, the model shows a damping error as shown in Figure 2a for the explicit method. Although Adam-Moulton's method is less accurate than the RK method, it still shows good stability analysis, as indicated in Figure 2e for values of ∆t ranging from 0-0.5. Thus, the Adam method is suitable for those who require frequent examination of the prevalence of diabetes to monitor minor changes. Note that the nature of Adam-Moulton's method depends on four previous values, as shown in Equations (40)-(45), to calculate the current value. Therefore, the values of the stability analysis of Adam-Moulton's method start after four steps.
Furthermore, the stability analysis of λ, γ, and more parameters in the diabetes model for the RK method was investigated. Figure 3 shows the analysis of the parameter γ (the rate of control of complications) which ranged between 0.04 and 1. it can be visualized that as the rate of complications increases, the number of diabetic patients with complications decreases. Hence, when gamma equals zero, meaning the complications are uncontrollable, death of diabetic people will occur due to the rise of diabetic patients with complications. However, if the rate of controlling the complications is slightly controlled, at a rate of 0.5, a substantial effect will be seen on the value of C. Moreover, for γ = 0.04, the value of C is more evened out. Ultimately, complications are always on the rise for diabetic patients since there is no permanent cure.
Consequently, Figure 4 shows the analysis for λ. The higher the probability of developing complications, the higher the number of diabetics with complications. It can be concluded from Figure 4, that the lower the probability of developing complications, diabetic patients without complications will increase. In developed countries, diabetes is an epidemic with most cases without complications. This is due to the low rate of developing complications and a high rate of controlling the complications through laws implemented by public health authorities [30]. Ultimately, our model shows consistent findings of the parametric analysis with Akinsola et al., who conducted a similar investigation [29].
In addition, the mortality of diabetic patients arises from the complications associated with the disease. In the DC model of diabetes, there are two parameters that are concerned with the mortality of diabetic patients with complications: the natural mortality rate (µ) and the mortality rate due to complications (δ). Figure 5 shows that the mortality rate due to complications increases, which corresponds to statistically more death cases of diabetic patients. Thus, the overall prevalence of diabetes in society decreases. Similarly, as the natural mortality rate increases in Figure 6, the value of C decreases.
Another parameter (ν) was analyzed using the RK method, as shown in Figure 7. This parameter is the rate at which patients with complications became severely disabled. People with diabetes-related disabilities [34] are at a higher risk of mortality as they are unable to control and manage their condition effectively. Therefore, this will increase the mortality rate due to diabetic complications, and eventually decrease the value of C. Care must be taken to reduce this rate by providing management programs to assist and support people in this category.
In Figure 8, the true value of C i is plotted in blue. However, it is not shown because it is very close to the values of the RK method. Thus, the graph distinctly illustrates that RK has the highest accuracy, followed by Adam-Moulton's method.
parametric analysis with Akinsola et al., who conducted a similar investigation [29].
In addition, the mortality of diabetic patients arises from the complications associated with the disease. In the DC model of diabetes, there are two parameters that are concerned with the mortality of diabetic patients with complications: the natural mortality rate ( ) and the mortality rate due to complications ( ). Figure 5 shows that the mortality rate due to complications increases, which corresponds to statistically more death cases of diabetic patients. Thus, the overall prevalence of diabetes in society decreases. Similarly, as the natural mortality rate increases in Figure 6, the value of C decreases.
Another parameter ( ) was analyzed using the RK method, as shown in Figure 7. This parameter is the rate at which patients with complications became severely disabled. People with diabetes-related disabilities [34] are at a higher risk of mortality as they are unable to control and manage their condition effectively. Therefore, this will increase the mortality rate due to diabetic complications, and eventually decrease the value of C. Care must be taken to reduce this rate by providing management programs to assist and support people in this category.
In Figure 8, the true value of Ci is plotted in blue. However, it is not shown because it is very close to the values of the RK method. Thus, the graph distinctly illustrates that RK has the highest accuracy, followed by Adam-Moulton's method.
Ultimately, for the RK method, the C values do not vary for different values of gamma, and this is due to the high accuracy of this method for the suggested ODE system. On the other hand, the Adam-Moulton method shows C values that are oscillating at γ equal to 0.5. This value was recommended by Boutayeb's paper [14]. However, it does not fit the model since it causes unstable results for the Adam-Moulton method. Hence, this paper proposes using γ to be equal to 0.08, which provides good stability for the Adam-Moulton method and the rest of the methods. For further analysis, the incidence of diabetes was considered in this paper. Three levels of incidence of diabetes are considered: low, medium, and high levels of incidence. Additionally, three different values of the parameters λ are studied for the various levels of incidence. By varying the parameter values, nine different events may be taken into consideration as can be seen in Table 6. For example, the scenario "high-high" denotes both a high rate of complications and a high incidence of diabetes. Ultimately, for the RK method, the C values do not vary for different values of gamma, and this is due to the high accuracy of this method for the suggested ODE system. On the other hand, the Adam-Moulton method shows C values that are oscillating at γ equal to 0.5. This value was recommended by Boutayeb's paper [14]. However, it does not fit the model since it causes unstable results for the Adam-Moulton method. Hence, this paper proposes using γ to be equal to 0.08, which provides good stability for the Adam-Moulton method and the rest of the methods.
For further analysis, the incidence of diabetes was considered in this paper. Three levels of incidence of diabetes are considered: low, medium, and high levels of incidence. Additionally, three different values of the parameters λ are studied for the various levels of incidence. By varying the parameter values, nine different events may be taken into consideration as can be seen in Table 6. For example, the scenario "high-high" denotes both a high rate of complications and a high incidence of diabetes. C is the number of diabetic patients with complications, N is the total number of people with diabetes, λ is the probability of developing a complication, and I is the incidence of diabetes.
From the above table, it can be noticed that there is an increase in both the longitudinal and latitudinal directions of the table. This indicates that the number of diabetics with complications and the total number of people with diabetes increase as the incidence of diabetes increases and as the rate of developing complications increases. This suggests that without effective health strategies and policies in place, the diabetes population may surge to a level beyond control, posing a threat to the resources of the country. In the situation of uncontrolled diabetes, the health authorities may not be able to serve or provide care for all the affected people, causing a significant threat to the socioeconomic welfare of a country.
The burden of diabetes can be significantly reduced if proper measures and schemes are employed at different levels to check the increase in diabetes incidence. First, by lowering the number of individuals who acquire prediabetes and individuals who get diabetes without complications, be it through diabetes education and health policies to promote good health and raise awareness, especially among the youth. There is a strong probability of reversing the effects of diabetes in its early stages. Secondly, by lowering the number of individuals with diabetes who have complications. This action has two main implications. Early detection of diabetes requires preventative measures, and after someone has been diagnosed with diabetes, care and precaution should be taken to prevent or, at the very least, delay the development of complications with the available resources [29].
Therefore, if the incidence is reduced, the number of pre-diabetics, diabetics, and diabetics with complications will also decrease. This goal can be realized through strategies that focus on raising awareness about a healthy diet, encouraging physical exercise, managing obesity, addressing blood pressure, reducing stress levels, quitting smoking and bad drinking habits, and regular screening.
We noted the ambiguity and dearth of data on the incidence of diabetes with complications and the worldwide status of diabetes complications rates. This gap in data on trends in complications has hindered the effective evaluation of the health situation of a country [35]. Future research on global trends in diabetes complications is necessary for better resource management and the identification of underlying factors and risks faced by the population.
The mathematical model illustration supports the diagnoses and recommendations of diabetes specialists and healthcare managers in general. Moreover, the results of our model are concurrent with the results from literation on similar models [19,22,29]. Thus, it provides health decision-makers with a framework for comparing the socioeconomic consequences of uncontrolled diabetes, and an effective strategy is a better investment in healthcare to reduce costs in the long run [19].

Conclusions
Diabetes and its complications pose an increasing threat to healthcare systems across the world. In this paper, a mathematical model of diabetes and the possible complications in a population was investigated and numerically analyzed. This paper added more numerical methods than were previously used. New methods such as the Heun's and the Adam-Moulton fourth-order methods were added. In addition, the step size (∆t) and its effect on the system were evaluated by computing the percentage error at every time-step size. The stability, numerical analysis, and simulations were determined for various stated values of the parameters for the analysis of the model. According to the results provided in the paper, the RK fourth-order method had the highest accuracy when compared with the other methods. Furthermore, by analyzing the parameters based on different variations, different situations can be generated, and subsequent strategies can be used in accordance with the available resource options. The optimal parameters' values were recommended to get the best results. From the graphs, if the incidence of diabetes is controlled, then the number of diabetics with or without complications also reduces.
In conclusion, controlling the rate of developing diabetes and monitoring complications in real life can be achieved through health education, healthy diet awareness, regular exercise, quitting smoking, and reducing other metabolic risks like obesity and hypertension. These actions will help in implementing optimal strategies to reduce the incidence of diabetes and the total number of diabetic patients. The numerical analysis of the model affirms the facts of rising diabetes incidence and prevalence around the world, according to IDF statistics. This emphasizes the significance of early diagnosis, monitoring, and treatment of diabetes mellitus to minimize complications, as well as provide diabetes patients with sufficient medical care. Special consideration should be given to policies and strategies that promote awareness of good health and the benefits of the prevention of diabetes rather than focusing primarily on curing the sick.
In the future work of this project, the partial differential equations can be investigated for more specific modeling of diabetes. Moreover, there should be more research with a special emphasis on finding the incidence of diabetes with and without complications. Research to find the incidence of diabetes showed that most international organizations do not report on the incidence of diabetes with complications. It is important to report this number for a better assessment of the health status of a country. To conclude, this paper can be improved by adding more mathematical methods to compare and make an affirmative conclusion in which method is highly recommended. Moreover, real-life diabetic data acquired from the public to be utilized for further testing of the model and provide attestation to our results. Lastly, if habits, health education and awareness, and diet are changed over the years, the results obtained from the model might differ significantly. the aforementioned points are considered as limitations in this paper.