A new modification of an iterative method based on inverse polynomial for solving Cauchy problems

Abstract

This paper introduces a new modification to an iterative method for solving Cauchy problems (IVPs) based on an inverse polynomial technique. The proposed method is proven to be consistent, stable, and convergent. We also demonstrate the consistency property of the One-stage scheme and the Two-stage scheme, as well as the stability property of the proposed method. To validate the accuracy of our approach, we conduct several numerical experiments and present the results graphically. Our findings show that the proposed method outperforms existing approaches in terms of accuracy and efficiency. Additionally, we discuss the implications of our results for future research in this field. Overall, this paper provides a valuable contribution to the numerical solution of IVPs and lays the groundwork for further exploration in this area.

Keywords:

• Approximation
• convergence
• differential equation
• initial value problem
• iterative method

1. Introduction

Differential equations play a crucial role in numerous fields of pure and applied sciences as they are employed to model a diverse range of real-world phenomena. Despite the availability of analytical methods for solving differential equations, many practical equations are too complex to have a closed-form solution. Even in cases where a solution formula exists, it may involve integrals that can only be approximated numerically. In such instances, numerical and iterative methods serve as an alternative approach to solving differential equations while satisfying specific initial conditions. Recently, several studies have been devoted to improving and modifying iterative numerical methods for efficiently, accurately, and reliably solving differential equations under diverse initial and boundary conditions. For instance, Ehiemua, Agbeboh, & Ataman (2020) proposed a new scheme based on the harmonic mean in 2020, while Olaniyan, Awe, & Akanbi (2021) established an improved third-order Runge-Kutta method utilizing the convex combination in 2021. In 2022, Yun (2022) presented a modification of the implicit Euler method for solving initial value problems. Additionally, Kadum & Abdul-Hassan (2023) conducted a recent study in 2023, introducing a novel method based on a quadrature integration formula. It is worth noting, however, that most of these studies did not consider the singularity problem in their experimental results. In this regard, inverse polynomial functions have been identified as a promising approach to address these challenging problems efficiently.

Inverse polynomial functions are used to develop a new method for solving differential equations that considers the equation’s particular character. This approach is based on the choice of the starting value of the function, which is selected to coincide with zero or the solution. This results in a technique that is smooth in the vicinity of singularities and can step over these points easily.

Recent studies have proposed several innovative approaches to solving initial value problems using numerical techniques. For example, in 2018, Fatma Rizaner and Ahmet Rizaner developed a method using unconstrained radial basis function networks to obtain approximate solutions for first-order initial value problems (Rizaner & Rizaner, 2018). In 2020, Kutniv et al introduced novel high-order one-step numerical techniques for initial value problems with a first-order singularity in the point $t=0$ (Kutniv, Datsko, Kunynets, & Włoch, 2020). Ramos and Rufai also developed a modified family of Falkner-type techniques for solving differential systems of second-order initial-value problems, which they analyzed and implemented using collocation and interpolation techniques (Ramos & Rufai, 2018, 2020).

The development of new techniques to solve higher-order boundary and initial value problems has been the subject of active research in the field of numerical analysis. Chebyshev-spectral method has been shown to be an efficient approach for solving such problems, as demonstrated in (Agarwal, Attary, Maghasedi, & Kumam, 2020), where the authors applied it to the several foundation problems. In addition, boundary shape functions methods have been used to solve nonlinear third-order three-point boundary value problems, as described in the work of Lin, Zhang, & Liu (2021). Another interesting approach is the symplectic integration of boundary value problems, which has been studied by several researchers, including McLachlan & Offen (2019). A comprehensive comparison of ODE solvers for biochemical problems has also been conducted to evaluate the performance of various numerical methods for solving such problems (Postawa, Szczygieł, & Kułażyński, 2020). Moreover, a new generalized operational matrix of integration based on Hermite wavelets has been introduced to solve nonlinear singular boundary value problems, as presented in the paper by Shiralashetti & Kumbinarasaiah (2019). Finally, boundary-domain integral method and homotopy analysis method have been used to solve systems of nonlinear boundary value problems in environmental engineering, as discussed in the work of Al-Jawary, Rahdi, & Ravnik (2020).

In this paper, we present a new modification of an iterative method based on an inverse polynomial technique for solving Cauchy problems. We prove the consistency, stability, and convergence of the proposed method and validate its accuracy through several numerical experiments. Our contribution adds to the growing body of research on numerical solution of differential equations and provides an innovative approach that can be utilized for a wide range of differential equations with diverse requirements.

As a result, and based on above motivations, we assumed in this study the regular Initial Value Problem (IVP) that takes the following form (1) $$u^{\prime }\left(t\right)=f\left(t,u\left(t\right)\right);\mathrm{ }u\left(t_0\right)=u_0.$$(1)

We consider here that the function $f\left(t,u\right)$ is globally Lipschitz continuous with respect to $u$ meaning that there exists a constant $L$ such that the Lipschitz condition $\left|f\left(t,u_1\right)-f\left(t,u_2\right)\right|\le L\left|u_1-u_2\right|$ is satisfied for all $t$ and $u_1,u_2,$ then the existence and uniqueness of solutions are guaranteed. Similarly, if $f(t,u)$ is continuously differentiable with respect to $u,$ then solutions exist and are unique in some neighborhood of the initial point (Ricardo, 2020).

One of the biggest challenges in solving differential equations is the singularities. Differential equations with singularities can be more challenging to solve and analyze than equations without singularities. In particular, singularities are points where the equation is not well-defined or behaves in an unusual way, such as points where the coefficient of the highest derivative of the dependent variable is zero or infinite. To deal effectively with differential equations with singularities, some particular methods can be taken into consideration. In these methods, the computations would lead to the inverse polynomial that takes the following form (2) $$u_{n+k}=\frac{u_n}{1+\sum _{j=1}^k\mathrm{\hspace{0.17em}}{b}_j{t}_n^j},$$(2) where the values of $b_j$ are real valued constants and $u_n$ is the numerical approximations to the solution at the point $t=t_n,$ and $k\ge 1$ is the step number and order of the methods.

The objectives of the study are to determine the numerical values of the parameters, $b_j$ for $k=1,2,$ computerized and implemented it on a microcomputer for numerical solution of some typical differential equations with singularities.

1.1. Determination of the coefficients of the methods

Setting $k=1$ in Equation (2)(2) $$u_{n+k}=\frac{u_n}{1+\sum _{j=1}^k\mathrm{\hspace{0.17em}}{b}_j{t}_n^j},$$(2) , we obtain a general one step formulae in the form (3) $$u_{n+1}=\frac{u_n}{1+b_1{t}_n}.$$(3)

Suppose that $\left|t\right|<1,$ then by the aid of binomial expansion theorem on the right hand side and ignoring terms of order higher than one, we have (4) $$u_{n+1}=u_n\left|1-b_1{t}_n\right|+O\left(t^2\right).$$(4)

The first order Taylor expansion of $u_{n+1}$ of point $\left(t_n,u_n\right)$ gives (5) $$u_{n+1}=u_n+h{u}_n^{\prime }+O\left(h^2\right).$$(5)

Substituting for $u_{n+1}$ in Equation (4)(4) $$u_{n+1}=u_n\left|1-b_1{t}_n\right|+O\left(t^2\right).$$(4) , we obtain (6) $$u_n+h{u}_n^{\prime }+O\left(h^2\right)=u_n-b_1{u}_n{t}_n+O\left(t^2\right).$$(6)

The coefficient $b_1$ is evaluated by imposing the condition that Equation (6)(6) $$u_n+h{u}_n^{\prime }+O\left(h^2\right)=u_n-b_1{u}_n{t}_n+O\left(t^2\right).$$(6) agrees term by term. Hence, we have (7) $$-b_1{u}_n{t}_n=h{u}_n^{\prime },$$(7) leading to (8) $$b_1=\frac{-h{u}_n^{\prime }}{u_n{t}_n}.$$(8)

Substituting this into Equation (3)(3) $$u_{n+1}=\frac{u_n}{1+b_1{t}_n}.$$(3) , we obtained a one-step integrator of first order in the form (9) $$u_{n+1}=\frac{u_n^2}{u_n-h{u}_n^{\prime }}.$$(9)

This incidentally coincides with inverse Euler method. The formula was used to approximate the solution of stiff and non-stiff ordinary differential equations numerically, and the results were found to be adequate.

Using the same argument above, by setting $k=2,$ in Equation (3)(3) $$u_{n+1}=\frac{u_n}{1+b_1{t}_n}.$$(3) , we obtain a second-order method as (10) $$u_{n+2}=\frac{u_n}{1+b_1{t}_n+b_2{t}_n^2}.$$(10)

Using again binomial and Taylor’s series expansion method (Abdul-Hassan, Ali, & Park, 2022; Ali & Pales, 2022) on Equation (10)(10) $$u_{n+2}=\frac{u_n}{1+b_1{t}_n+b_2{t}_n^2}.$$(10) , we obtain (11) $$u_n+h{u}_n^{\prime }+\frac{h^2}{2}{u}_n^{″}+O\left(h^3\right)=u_n-u_n{b}_1{t}_n-u_n{b}_2{t}_n^2+u_n{b}_1^2{t}_n^2{u}_n^2+O\left(t^3\right).$$(11)

Using term-by-term equality principle in Equation (11)(11) $$u_n+h{u}_n^{\prime }+\frac{h^2}{2}{u}_n^{″}+O\left(h^3\right)=u_n-u_n{b}_1{t}_n-u_n{b}_2{t}_n^2+u_n{b}_1^2{t}_n^2{u}_n^2+O\left(t^3\right).$$(11) . We obtain the set of non-linear algebraic equations (12) $$\left.\begin{array}{c}-b_1{u}_n{t}_n=h{u}_n^{\prime },\\ -t_n^2\left(u_n{b}_2-u_n{b}_1^2\right)=\frac{h^2{u}_n^{″}}{2}\end{array}\right\}.$$(12)

Solving these set of algebraic equations, we get (13) $$b_1=\frac{-h{u}_n^{\prime }}{u_n{t}_n},\mathrm{ and}\mathrm{ }{b}_2=\frac{h^2}{2u_n^2{t}_n^2}\left[2{\left(u_n^{\prime }\right)}^2-u_n{u}_n^{″}\right].$$(13)

Substituting for $b_1$ and $b_2$ in Equation (10)(10) $$u_{n+2}=\frac{u_n}{1+b_1{t}_n+b_2{t}_n^2}.$$(10) , we obtain a two-stage method of second-order as the following (14) $$u_{n+2}=\frac{2u_n^3}{2u_n^2-2h{u}_n{u}_n^{\prime }+h^2\left(2{\left(u_n^{\prime }\right)}^2-u_n{u}_n^{″}\right)}.$$(14)

In the next sections, we study the consistency and stability properties of the above formulas.

2. Properties of the formula

2.1. Consistency property

2.1.1. One-stage scheme

Subtracting $u_n$ from both sides of Equation (9)(9) $$u_{n+1}=\frac{u_n^2}{u_n-h{u}_n^{\prime }}.$$(9) , we get (15) $$u_{n+1}-u_n=\frac{h{u}_n{u}_n^{\prime }}{u_n-h{u}_n^{\prime }}.$$(15)

Dividing both sides by $h,$ one could get (16) $$\frac{u_{n+1}-u_n}{h}=\frac{u_n{u}_n^{\prime }}{u_n-h{u}_n^{\prime }}.$$(16)

Taking the limit, as $h$ tends to zero, on both sides of Equation (16)(16) $$\frac{u_{n+1}-u_n}{h}=\frac{u_n{u}_n^{\prime }}{u_n-h{u}_n^{\prime }}.$$(16) , we have (17) $$\underset{h\to 0}{\lim }\left(\frac{u_{n+1}-u_n}{h}\right)\to u_n^{\prime }=f\left(t_n,u_n\right).$$(17)

Suggesting that the one-stage integrator in Equation (9)(9) $$u_{n+1}=\frac{u_n^2}{u_n-h{u}_n^{\prime }}.$$(9) is consistent.

2.1.2. Two-stage scheme

Similarly, subtracting $u_n$ from both sides of Equation (14)(14) $$u_{n+2}=\frac{2u_n^3}{2u_n^2-2h{u}_n{u}_n^{\prime }+h^2\left(2{\left(u_n^{\prime }\right)}^2-u_n{u}_n^{″}\right)}.$$(14) and dividing by $2h,$ we get (18) $$\frac{u_{n+2}-u_n}{2h}=\frac{u_n^2{u}_n^{\prime }-\frac{h}{2}{u}_n^{\prime }\left(2{\left(u_n^{\prime }\right)}^2-u_n{u}_n^{″}\right)}{u_n^2-h{u}_n{u}_n^{\prime }+\frac{h^2}{2}\left(2{\left(u_n^{\prime }\right)}^2-u_n{u}_n^{″}\right)}.$$(18)

Taking limit on both sides of (19)(19) $$\underset{h\to 0}{\lim }\left(\frac{u_{n+2}-u_n}{2h}\right)\to u_n^{\prime }.$$(19) as $h\to 0,$ we have (19) $$\underset{h\to 0}{\lim }\left(\frac{u_{n+2}-u_n}{2h}\right)\to u_n^{\prime }.$$(19)

Implying that the two-stage scheme in Equation (14)(14) $$u_{n+2}=\frac{2u_n^3}{2u_n^2-2h{u}_n{u}_n^{\prime }+h^2\left(2{\left(u_n^{\prime }\right)}^2-u_n{u}_n^{″}\right)}.$$(14) is consistent.

2.2. Stability property

Applying the numerical integrator of Equation (14)(14) $$u_{n+2}=\frac{2u_n^3}{2u_n^2-2h{u}_n{u}_n^{\prime }+h^2\left(2{\left(u_n^{\prime }\right)}^2-u_n{u}_n^{″}\right)}.$$(14) to the following stability test equation (20) $$\begin{array}{ll}& u^{\prime }=\lambda u;u\left(t_0\right)=u_0,\mathrm{ }Re\left(\lambda \right)<0\\ & a\le t\le b\end{array},$$(20) we obtain a finite difference equation (21) $$u_{n+1}=\left[\frac{1}{1-\lambda h}\right]u_n,$$(21) with $\frac{1}{1-\lambda h}$ as the stability function.

For the convergence of the scheme, we have (22) $$\left|\mu \left(z\right)\right|<1,\mathrm{ }z=\lambda h.$$(22)

Analysis of Equation (22)(22) $$\left|\mu \left(z\right)\right|<1,\mathrm{ }z=\lambda h.$$(22) gives $z<0,$ showing that the corresponding interval of absolute stability of the scheme is $\left(-\infty ,0\right).$

Similarly, applying the formula of Equation (14)(14) $$u_{n+2}=\frac{2u_n^3}{2u_n^2-2h{u}_n{u}_n^{\prime }+h^2\left(2{\left(u_n^{\prime }\right)}^2-u_n{u}_n^{″}\right)}.$$(14) to Equation (20)(20) $$\begin{array}{ll}& u^{\prime }=\lambda u;u\left(t_0\right)=u_0,\mathrm{ }Re\left(\lambda \right)<0\\ & a\le t\le b\end{array},$$(20) yields (23) $$u_{n+2}=\left[\frac{1}{1-\lambda h+\frac{{\lambda }^2\mathrm{ }{h}^2}{2}}\right]u_n,$$(23) which leads to a convergent solution, whenever the stability function $\mu \left(\lambda h\right)$ which is given by (24) $$\mu \left(\lambda h\right)=\frac{1}{1-\lambda h+\frac{{\lambda }^2{h}^2}{2}}$$(24) satisfies the following inequality: (25) $$\left|\mu \left(\lambda h\right)\right|<1.$$(25)

Setting $z=\lambda h,$ then the absolute stability interval of the scheme is seen to lie in the Region defined by set (26) $$\beta =z:\left|\mu \left(\lambda h\right)\right|<1.$$(26)

That is, ${\beta }_1=z:z<0$ and ${\beta }_2=z:z<2.$ Showing that the two-stage scheme is weakly stable and not suitable for stiff equations. Since the scheme is consistent and weakly stable then it is convergent whenever $\lambda h\in \beta .$

3. Numerical experiments

To evaluate the efficiency and accuracy of our proposed method, we implemented the numerical formula using the MATLAB programming language and tested it with three examples. Examples 1 and 2 demonstrate the performance of the method on a well-known low-order differential equation with a nonsingular solution, while Example 3 demonstrates the performance of the method on a differential equation with a point of singularity. To further validate the effectiveness of the method, we compared it with three well-known methods, Heun’s Method (Ascher & Petzold, 1998), Ralston’s Method (Ralston, 1962), and the Runge-Kutta Method (RK2 Method) (Hatun & Vatansever, 2016), which have the same order of convergence. We used various step sizes to ensure the validity of the tests.

3.1. Example 1

Consider the non-singular equation, given by $$u^{\prime }=u-t^2+1;\mathrm{ }u\left(0\right)=0.5.$$

The exact solution is given by $$u\left(t\right)=(1+t{)}^2-\frac{e^t}{2}.$$

The results of the comparison with the related methods are presented in Table 1. Additionally, the absolute error is visualized in Figure 1.

Figure 1. The absolute error in Example 1 of the proposed method and other related methods.

Table 1. Comparison in Example 1 of the proposed method, related methods, and the exact solution.

$t$	Proposed method	Heun’s method	Ralston’s method	RK2 method	Exact solution
0	0.5	0.5	0.5	0.5	0.5
0.04	0.561184000000000	0.561168000000000	0.561176000000000	0.561178666666667	0.561194612903806
0.08	0.624735027200000	0.624702374400000	0.624718700800000	0.624724142933333	0.624756466162521
0.12	0.690619096309760	0.690569111275520	0.690594103792640	0.690602434631680	0.690651574210312
0.16	0.758800835439198	0.758732811015561	0.758766823227380	0.758778160631319	0.758844564504095
0.2	0.829243429525118	0.829156629704996	0.829200029615057	0.829214496251744	0.829298620919915
0.24	0.901908561449742	0.901802220196960	0.901855390823351	0.901873114365482	0.901975424839298
0.28	0.976756350756892	0.976629670780996	0.976693010768944	0.976714124098260	0.976835093831282
0.32	1.05374528986777	1.05359744134886	1.05367136560832	1.05369600702814	1.05383611783202
0.36	1.13283217769438	1.13266229695589	1.13274723732514	1.13277555078155	1.13293529271983
0.4	1.21397205054431	1.21377923867169	1.21387564460800	1.21390777992010	1.21408765117936
0.44	1.29711811020652	1.29690143160950	1.29700977090801	1.29704588400751	1.29724639074433
0.48	1.38222164910294	1.38198013001917	1.38210088956106	1.38214114274168	1.38236279890355
0.52	1.46923197238634	1.46896459932395	1.46909828585515	1.46914284803221	1.46938617515056
0.56	1.55809631685971	1.55780203497637	1.55794917591804	1.55799822289859	1.55826374985195
0.6	1.64875976658758	1.64843747800340	1.64859862229549	1.64865233705952	1.64894059980475
0.64	1.74116516506435	1.74081372710594	1.74098944608515	1.74104801907822	1.74135956034752
0.68	1.83525302379898	1.83487124717186	1.83506213548542	1.83512576492327	1.83546113388478
0.72	1.93096142716998	1.93054807405648	1.93075475061323	1.93082364279881	1.93118339467806
0.76	2.02822593339851	2.02777971547798	2.02800282443825	2.02807719409167	2.02846188975159
0.8	2.12697947148117	2.12649904786948	2.12673925967533	2.12681933027728	2.12722953575377
0.84	2.22715223391761	2.22663620902256	2.22689422147008	2.22698022561926	2.22741651160946
0.88	2.32867156506144	2.32811848635068	2.32839502570606	2.32848720549119	2.32895014679140
0.92	2.43146184491595	2.43087020059378	2.43116602275487	2.43126463014189	2.43175480503185
0.96	2.53544436818852	2.53481258477801	2.53512847648327	2.53523377371835	2.53575176328844
1	2.64053721841061	2.63986365823695	2.64020043832378	2.64031269835273	2.64085908577048

3.2. Example 2

Consider the following initial value problem $$u^{\prime }=e^{t-u};\mathrm{ }u\left(0\right)=1,$$ providing that the exact solution is $$u\left(t\right)=ln\left(e^t+e-1\right).$$

The results of the comparison with the related methods are presented in Table 2. Additionally, the absolute error is visualized in Figure 2.

Figure 2. The absolute error in Example 2 of the proposed method and other related methods.

Table 2. Comparison in Example 2 of the proposed method, related methods, and the exact solution.

$t$	Proposed method	Heun’s method	Ralston’s method	RK2 method	Exact solution
0	1	1	1	1	1
0.04	1.01490239391214	1.01482715726501	1.01490298812788	1.01490278977730	1.01490185851801
0.08	1.03018067795801	1.03002943061651	1.03018185449769	1.03018146177012	1.03017958030835
0.12	1.04583841779603	1.04561043590589	1.04584016385853	1.04583958103103	1.04583673137268
0.16	1.06187892229230	1.06157353435455	1.06188122419522	1.06188045583891	1.06187662103953
0.2	1.07830523368655	1.07792182259499	1.07830807689868	1.07830712786790	1.07830229213246
0.24	1.09512011843491	1.09465812337622	1.09512348761310	1.09512236303308	1.09511651181413
0.28	1.11232605878669	1.11178497699116	1.11232993781694	1.11232864307005	1.11232176316363
0.32	1.12992524514884	1.12930463348037	1.12992961719127	1.12992815790202	1.12992023754031
0.36	1.14791956928694	1.14721904566202	1.14792441682460	1.14792279884331	1.14791382778352
0.4	1.16631061840743	1.16552986303358	1.16631592329847	1.16631415268384	1.16630412229232
0.44	1.18509967015974	1.18423842658516	1.18510541369275	1.18510349669338	1.18509240002431
0.48	1.20428768859184	1.20334576455906	1.20429385154395	1.20429179457894	1.20427962644680
0.52	1.22387532108615	1.22285258918376	1.22388188378351	1.22387969342224	1.22386645046740
0.56	1.24386289629646	1.24275929440448	1.24386983867673	1.24386752161802	1.24385320236481
0.6	1.26425042309995	1.26306595462543	1.26425772477633	1.26425528782705	1.26423989273362
0.64	1.28503759057119	1.28377232447255	1.28504523089765	1.28504268095092	1.28502621245045
0.68	1.30622376897854	1.30487783957806	1.30623172711563	1.30622907112887	1.30621153366164
0.72	1.32780801179587	1.32638161838174	1.32781626677693	1.32781351174978	1.32779491178561
0.76	1.34978905871636	1.34828246493655	1.34979758951338	1.34979474246566	1.34977508851664
0.8	1.37216533964760	1.37057887269981	1.37217412523657	1.37217119318652	1.37215049580954
0.84	1.39493497966160	1.39326902928446	1.39494399908671	1.39494098902954	1.39491926081876
0.88	1.41809580486645	1.41635082213864	1.41810503730288	1.41810195618983	1.41807921175881
0.92	1.44164534916108	1.43982184511623	1.44165477397590	1.44165162869390	1.44162788464772
0.96	1.46558086182892	1.46367940589537	1.46559045863979	1.46558725599189	1.46556253088918
1	1.48989931592156	1.48792053419692	1.48990906465286	1.48990581133951	1.48988012564475

3.3. Example 3

Consider the following initial value problem $$u^{\prime }=\frac{u\mathrm{ }\log \left(u\right)}{1-t};\mathrm{ }u\left(0\right)=e^{0.2},$$ providing that the exact solution is $$u\left(t\right)=e^{\frac{0.2}{1-t}},$$ with singularity at $t=1.$

The results of the comparison with the related methods are presented in Table 3. Additionally, the absolute error is visualized in Figure 3 (from $t=0$ to $t=0.48$) and in Figure 4 (from $t=0\mathrm{ }$ to $\mathrm{ }t=0.84$).

Figure 3. The absolute error in Example 3 of the proposed method and other related methods ($t=0\mathrm{ }to\mathrm{ }t=0.48$).

Figure 4. The absolute error in Example 3 of the proposed method and other related methods ($t=0\mathrm{ }to\mathrm{ }t=0.84$).

Table 3. Comparison in Example 3 of the proposed method, related methods, and the exact solution.

$t$	Proposed method	Heun’s Method	Ralston’s method	RK2 method	Exact solution
0	$e^{0.2}$	$e^{0.2}$	$e^{0.2}$	$e^{0.2}$	$e^{0.2}$
0.04	1.23162264000921	1.23161308644960	1.23161781411736	1.23161622746386	1.23162364234705
0.08	1.24282802421235	1.24280653210221	1.24281716501400	1.24281359592009	1.24283033191213
0.12	1.25516812370438	1.25513163624676	1.25514968288116	1.25514362416900	1.25517214093937
0.16	1.26882375585456	1.26876831582900	1.26879572788691	1.26878652311555	1.26883002802581
0.2	1.28401614396063	1.28393657116747	1.28397590191972	1.28396269196626	1.28402541668774
0.24	1.30101882127246	1.30090824985149	1.30096288108010	1.30094452753854	1.30103212886034
0.28	1.32017398754713	1.32002319004898	1.32009766326784	1.32007263658142	1.32019278843412
0.32	1.34191538091227	1.34171176056579	1.34181227030819	1.34177848305191	1.34194176988288
0.36	1.36680089182732	1.36652694562098	1.36666209086979	1.36661664359580	1.36683794117380
0.4	1.39556010096327	1.39519101505220	1.39537297331491	1.39531175710960	1.39561242508609
0.44	1.42916531142209	1.42866508873993	1.42891150390721	1.42882856070453	1.42924003241798
0.48	1.46894073167291	1.46825574686909	1.46859286578793	1.46847932439096	1.46904919384902
0.52	1.51673584760925	1.51578370555336	1.51625178914067	1.51609402733599	1.51689679638821
0.56	1.57521132277439	1.57386073041423	1.57452380593205	1.57430013317472	1.57545710339032
0.6	1.64833188580287	1.64636435357128	1.64732871763507	1.64700306992639	1.64872127070013
0.64	1.74226230743981	1.73929462769733	1.74074616880987	1.74025536332157	1.74290899863346
0.68	1.86710425861103	1.86241936828633	1.86470464613910	1.86393061998811	1.86824595743222
0.72	2.04054218373681	2.03268423601039	2.03650352177981	2.03520701455561	2.04272707026614
0.76	2.29631086834141	2.28199453506674	2.28891798294771	2.28656043348219	2.30097589089283
0.8	2.70665442159724	2.67733590010588	2.69140989038687	2.68659585672887	2.71828182845905
0.84	3.45380795556122	3.38236145498047	3.41625657675448	3.40458117696516	3.49034295746185
0.88	5.12702310648509	4.89723834782606	5.00398744172861	4.96677151210114	5.29449005047005
0.92	10.6214628415599	9.41525489116971	9.95129370042096	9.75989689757794	12.1824939607036
0.96	53.3059166397085	35.7396518978688	42.6045493713402	40.0029551487507	148.413159102582
1	Diverges	Diverges	Diverges	Diverges	0

4. Results and discussion

Tables 1–3 present the results of applying the proposed method to Example 1, Example 2, and Example 3, respectively, with a step size of $h=0.04.$ The tables show the evaluation point ($t$) in the first column, followed by the corresponding approximate solutions obtained by the proposed method, Heun’s method, Ralston’s method, and the Runge-Kutta method in the second to fifth columns. The last column of each table displays the exact solution of the initial value problem, computed using the given exact function in each example. On the other hand, we calculate the absolute error is using the following formula: $$\mathit{Absolute}\mathrm{ }\mathit{error}=\left|\mathit{Exact}\mathrm{ }\mathit{solution}\mathrm{ }\text{-}\mathrm{ }\mathit{Approximated}\mathrm{ }\mathit{solution}\right|.$$

The formula is evaluated for each point in the solution and the resulting values are presented and illustrated in Figures 1, 2 for Example 1 and Example 2, respectively, while we had to avoid visualizing Example 3 on the entire interval $[0,1]$ due to the singularity at 1 as well as the exponential divergence from the solution when we approach to the point of singularity. Therefore, we illustrated in Figures 3, 4 the absolute error of Example 3 (from $t=0$ to $t=0.48$) to see the performance of the proposed method clearly, and (from $t=0\mathrm{ }$ to $\mathrm{ }t=0.84$) to see the overall results, respectively.

The numerical simulations and comparisons presented in this paper were conducted using MATLAB (R2018a) on a personal computer equipped with Windows 10 Pro, a 10th Generation Intel(R) Core (TM) i7-10600H @ 2.50 GHz processor, and 8.0 GB RAM storage. The superiority of the proposed method over the other methods was evident in the three examples. The results show that the proposed method is highly accurate, and its performance is excellent compared to other methods, as shown in Tables 1–3. Furthermore, the accuracy of the proposed method was demonstrated through absolute error analysis, which confirmed its efficiency. In Examples 1–3, we tested the accuracy of the proposed method by taking a step size of $h=0.04.$ The results show that the proposed method yields highly accurate approximations, and its accuracy improves as the step size is reduced.

5. Conclusions

In this study, we introduced a new modification of an iterative method based on inverse polynomial techniques for the numerical solution of Cauchy problems. Our proposed two-step second-order approach was demonstrated to be consistent, stable, and effective in solving both singular and non-singular ordinary differential equations through various computational experiments. The numerical results from our proposed method confirmed its accuracy and effectiveness, highlighting its versatility and ability to solve problems arising from various fields, including electrical networks, inflation-affected economies, and chemical reactions.

However, it is worth noting that the accuracy of the proposed method can be further enhanced by using Richardson error control, which we recommend as an avenue for future research. Overall, this paper contributes to the advancement of efficient numerical methods for the solution of ordinary differential equations and has the potential to facilitate scientific research in numerous fields.

Disclosure statement

All authors have declared they do not have any competing interests.

Data availability statement

No data were used to support this study.

Additional information

Funding

This research did not receive any external funds.