Exploring a Novel Multi-Stage Differential Transform Method Coupled with Adomian Polynomials for Solving Implicit Nonlinear ODEs with Analytical Solutions

Abstract

In engineering, physics, and other fields, implicit ordinary differential equations are essential to simulate complex systems. However, because of their intrinsic nonlinearity and difficulty separating higher-order derivatives, implicit ordinary differential equations pose substantial challenges. When applied to these types of equations, traditional numerical methods frequently have problems with convergence or require a significant amount of computing power. In this work, we present the multi-stage differential transform method, a novel semi-analytical approach for effectively solving first- and second-order implicit ordinary differential systems, in conjunction with Adomian polynomials. The main contribution of this method is that it simplifies the solution procedure and lowers processing costs by enabling the differential transform method to be applied directly to implicit systems without transforming them into explicit or quasi-linear forms. We obtain straightforward and effective algorithms that build solutions incrementally utilizing the characteristics of Adomian polynomials, providing benefits in theory and practice. By solving several implicit ODE systems that are difficult for traditional software programs such as Maple 2024, Mathematica 14, or Matlab 24.1, we validate our approach. The multi-stage differential transform method’s contribution includes expanded convergence intervals for numerical results, more accurate approximate solutions for wider domains, and the efficient calculation of exact solutions as a convergent power series. Because of its ease of implementation in educational computational tools and substantial advantages in terms of simplicity and efficiency, our method is suitable for researchers and practitioners working with complex implicit differential equations.

1. Introduction

Differential equations are essential for modeling a broad range of phenomena in physics, engineering, biology, and economics, among other disciplines. To forecast future behavior based on existing conditions, these equations offer a mathematical tool for characterizing dynamical systems and processes. For explicit ODE systems, numerous solution methods are available, including the Runge–Kutta method [1,2], the power series method [3], the Adomian decomposition method [4,5,6,7], the homotopy perturbation method [8], and the differential transform method (DTM) [9,10,11]. In addition, many software packages are capable of solving these explicit systems. Because of their capacity to represent complex systems in which the relationship between variables is not explicitly defined, implicit ordinary differential equations (IODEs) have attracted much interest among the various kinds of differential equations. In these systems, the higher derivative of the unknown cannot, in general, be easily isolated due to nonlinearities. Although these systems can offer greater flexibility in modeling complex problems, they can be more challenging to solve and require more advanced numerical methods. Packages such as Mathematica 14, Matlab 24.1, and Maple 2024 are not able to solve these systems in many situations, and examples of such cases are provided. Furthermore, applying implicit numerical integration methods to such systems will require the solution of nonlinear algebraic systems at every step. This will involve a large computational effort, especially for large-dimension systems.

Implicit differential equations form a major topic of study in the science of dynamical systems since they emerge in engineering situations. IODEs involve several major difficulties compared to explicit ODEs. First, the existence and uniqueness of solutions can be challenging. Second, the numerical integration of IODEs is more complex; conventional approaches sometimes depend on converting these equations into explicit forms, which can be computationally expensive. In addition, IODEs often have nonlinearities that could cause numerical method convergence problems, therefore complicating the process of finding precise responses. These elements call for the development of specific methods for successfully treating IODEs in different uses [12].

Furthermore, PDEs and ODEs occur in many technological and scientific domains. Even though many of these equations are PDEs, they are frequently converted to ODEs using traveling waves or other standard reductions [13]. Some nonlinearities in the PDEs lead to implicit ordinary differential equations.

Forces that depend on accelerations produce IODEs in mechanical systems, especially those involving dry friction, which are extremely complex. Because of these dependencies, implicit differential equations cannot be solved for higher derivatives, creating significant mathematical difficulties about the existence and uniqueness of solutions. Dry friction in sliding and rolling mechanical systems can behave non-deterministically, producing multiple solutions with the same starting conditions—a phenomenon associated with the Painlevé paradoxes. Furthermore, dry friction can appear as viscous resistance proportional to the square of velocity due to implicit equations in these systems, which can also produce stagnation zones where motion stops. Since implicit equations have considerable effects on system dynamics and might produce unexpected behaviors, they must be carefully considered while studying mechanical systems. To effectively forecast and control the motion of complex mechanical systems, a deep comprehension of these equations is essential [14].

Through index analysis and computational differentiation, differential algebraic equations (DAEs) can frequently be converted into equivalent IODEs under specific circumstances, which is particularly interesting given that DAEs have applications in fields such as control theory, fluid dynamics, circuit analysis, and mechanical and chemical engineering [12,15,16,17]. Among these diverse engineering applications, the field of electrical engineering provides particularly compelling examples of IODEs. The analysis of nonlinear RLC circuits is one of the most important uses of these equations since it offers a framework for comprehending the complex dynamics that result from the interaction of resistance, inductance, and capacitance. In these circuits, the system’s behavior is frequently characterized by equations that are not explicitly solvable for the higher derivatives, hence implicit formulations are required [16].

Differential systems of the type $\mathbf{F}\left({\mathbf{x}}^{\left(m\right)}\left(t\right),\dots ,{\mathbf{x}}^{\prime }\left(t\right),\mathbf{x}\left(t\right),t\right)=\mathbf{0}$, where the Jacobian of $\mathbf{F}$ with respect to the mth derivative ${\mathbf{x}}^{\left(m\right)}$ is nonsingular, are known as implicit ordinary differential systems. These equations require more complex methods for analysis and solution than explicit ones, where the mth derivative can be isolated. Implicit relationships complicate the solution procedure because they may not be immediately applicable to standard methods [16,17]. Consequently, several numerical methods for finding approximate solutions have been developed through research on IODEs.

By parameterizing IODEs, researchers can describe the equations in a way that shows the connections among variables without having to rewrite the system in an explicit form. Through this parametrization, researchers can obtain general solutions by converting implicit ODEs into standard systems of ordinary differential equations. These results can then be used to solve related differential algebraic problems [18]. While this parametrization approach provides a theoretical framework, the practical implementation of solutions requires specific numerical methods. A variety of techniques have been proposed to solve explicit ODEs. Runge–Kutta methods, for example, are frequently employed due to their efficiency in resolving ordinary differential equations. In contrast, they encounter difficulties with implicit equations as a result of the need for iterative solutions, which can result in intensified computational costs and convergence issues [13,19]. Differential equations can be transformed into algebraic equations using finite difference methods. Despite their effectiveness, these techniques may be computationally costly since they need exact discretization and have stability problems [13,19]. In addition, Picard iteration and incremental harmonic balance methods have been implemented; nevertheless, they can fail and commonly face convergence issues [13,19].

From theoretical convergence frameworks to useful numerical techniques and visualizations, the study of implicit ODEs uses a variety of approaches. To guarantee the presence of unique solutions, for instance, consecutive approximations in Banach spaces offer an organized method using operator techniques [20]. By using the Vessiot distribution, the geometric theory reduces implicit ODEs to ordinary dynamical systems analysis by transforming them into vector field distributions on manifolds [21]. Techniques that emphasize convergence conditions, such as the backward Euler method, are used to numerically solve initial and boundary value problems [22]. The Runge–Kutta and multi-step methods for implicit equations are also modified to maintain the same level of accuracy as explicit ODEs [23]. The mathematical basis of quadrature methods is further expanded by introducing an existence theorem and providing examples [24].

Given the importance of IODEs and the aforementioned numerical issues, we introduce a new semi-analytical method called the multi-stage differential transform method (MsDTM), coupled with Adomian polynomials. This approach presents a systematic means of obtaining implicit ODE solutions by constructing solutions incrementally. The following summarizes the MsDTM. The MsDTM subdivides the solution interval into a set number of sub-intervals for its operation. In the first sub-interval, the initial condition enables the first application of the DTM to the differential equation, producing the solution for this interval. Next, the algorithm processes the second sub-interval. Again, the DTM applies to the differential equation using the initial condition now derived from the first sub-interval’s solution. This step ensures continuity and consistency of the solution between sub-intervals. For each subsequent sub-interval, the same procedure is repeated. For each case, the prior sub-interval’s solution updates the next sub-interval’s initial condition. As we follow the algorithm, this iterative process yields a progressive solution. When the algorithm reaches the final step, it finds a solution for the whole interval. Finally, each sub-interval’s contribution ensures the method accurately approximates the solution across the entire domain. Readers can find more details in references [25,26,27]. The ease of use of MsDTM in computational tools such as Maple, Matlab, or Mathematica is one of the main advantages since it allows access to both scholars and practitioners. The MsDTM, combined with Adomian polynomials, seeks to provide precise solutions to IODEs by utilizing the advantages of both the differential transform method and Adomian polynomials.

Then, we propose a novel and efficient technique to solve implicit mth order ($m=1,m=2$) systems of ODEs of the form:

$$\mathbf{F}\left({\mathbf{x}}^{\left(m\right)}\left(t\right),\dots ,{\mathbf{x}}^{\prime }\left(t\right),\mathbf{x}\left(t\right),t\right)=\mathbf{0},t\in [0,T],$$

(1)

where $\mathbf{x}\left(t\right)\in {\mathbb{R}}^n$ is the sought solution and ${\mathbf{x}}^{\left(m\right)}\left(t\right)$ denotes its mth derivative and $\mathbf{x}:{\left({\mathbb{R}}^n\right)}^{(m+1)}\times [0,T]⟶{\mathbb{R}}^n$. This system is subject to the following initial conditions:

$${\mathbf{x}}^{\left(i\right)}\left(0\right)={\mathbf{x}}_i,0\le i\le m-1,$$

(2)

where ${\mathbf{x}}_i\in {\mathbb{R}}^n$ are given vectors. Throughout this paper, we assume the function $\mathbf{F}({\mathbf{y}}_m,\dots ,{\mathbf{y}}_1,{\mathbf{y}}_0,t)$ to be analytical with respect to its arguments, and the Jacobian $J:={\mathbf{F}}_{{\mathbf{y}}_m}^{\prime }({\mathbf{y}}_m,\dots ,{\mathbf{y}}_1,{\mathbf{y}}_0,t)\in {\mathbb{R}}^{n\times n}$ to be nonsingular. Furthermore, the present work focuses on the numerical method for solving IODEs. Due to the complexity of the IODEs discussed in this paper, we assume the solutions of these systems to be analytical. Some sub-classes of (1) are of the following explicit form:

$${\mathbf{x}}^{\left(m\right)}\left(t\right)=\mathbf{G}\left({\mathbf{x}}^{(m-1)}\left(t\right),\dots ,{\mathbf{x}}^{\prime }\left(t\right),\mathbf{x}\left(t\right),t\right),t\in [0,T],$$

(3)

and of the quasi-linear implicit form:

$$\mathbf{A}\left({\mathbf{x}}^{(m-1)}\left(t\right),\dots ,{\mathbf{x}}^{\prime }\left(t\right),\mathbf{x}\left(t\right),t\right){\mathbf{x}}^{\left(m\right)}\left(t\right)=\mathbf{G}\left({\mathbf{x}}^{(m-1)}\left(t\right),\dots ,{\mathbf{x}}^{\prime }\left(t\right),\mathbf{x}\left(t\right),t\right),t\in [0,T],$$

(4)

where $\mathbf{A}:=\mathbf{A}\left({\mathbf{x}}^{(m-1)}\left(t\right),\dots ,{\mathbf{x}}^{\prime }\left(t\right),\mathbf{x}\left(t\right),t\right)\in {\mathbb{R}}^{n\times n}$ is a nonsingular matrix.

The technique we propose combines the DTM with Adomian polynomials [28,29,30,31,32,33] to create a powerful solution method. The key idea is to apply the DTM directly to the system (1) without converting it to the semi-explicit form (3) or the quasi-linear form (4). By leveraging an important property of Adomian polynomials, we have derived a simple and efficient algorithm for this new technique, significantly enhancing its simplicity and effectiveness.

Our approach can yield the exact solution as a convergent power series, provided that all computations are performed with exact precision. To extend the convergence interval of the DTM solution, we have developed a multi-stage DTM algorithm (MsDTM). Both, the DTM and the MsDTM algorithms can be easily implemented using software packages such as Maple, Mathematica, or Matlab.

To demonstrate the effectiveness and accuracy of our technique, we solved several implicit systems of ODEs using the new method. It is important to note that all the numerical examples presented in this paper are not solvable by Maple 2024, numerically and symbolically. Furthermore, the Runge–Kutta method is also inapplicable to these examples. The numerical results show that our technique has successfully solved all given numerical examples. Specifically, the standard DTM provides the exact solution as a convergent power series, while the MsDTM produces accurate approximate solutions over larger intervals. To the best of our knowledge, this paper is the first to explore the numerical solution of implicit ordinary differential systems.

This paper’s main contribution is a proposed MsDTM that solves IODEs without needing transformation to explicit or quasi-linear forms, thus reducing computational costs. Adomian polynomials, when integrated, expand nonlinear functions into power series to solve difficult problems. The method, unlike traditional Runge–Kutta methods, handles first- and second-order implicit systems and generates accurate solutions expressed as convergent power series. With high accuracy, the multi-stage approach extends convergence intervals.

We structure the paper as follows: In Section 2, we review Adomian polynomials and the DTM. Section 3 presents two theorems that form the basis of the new algorithms for solving the implicit initial-value problems (1) and (2). In Section 4, we solve several first and second-order implicit systems of ODEs with different types of nonlinearities to illustrate the accuracy of the new method. In Section 5, we focus on the numerical procedure and examples. Finally, in Section 6, we present our concluding remarks.

2. Adomian Polynomials and the DTM

In this section, we give a brief review of Adomian polynomials [28,29,30,31,32,33] and the differential transform method. The differential transforms of the functions $\mathbf{x}\left(t\right)$ and $\mathbf{y}\left(t\right)$ are defined by the following:

$$x_k={\displaystyle \frac{1}{k!}}{\left[{\displaystyle \frac{d^k\mathbf{x}\left(t\right)}{d{t}^k}}\right]}_{t=0},y_k={\displaystyle \frac{1}{k!}}{\left[{\displaystyle \frac{d^k\mathbf{y}\left(t\right)}{d{t}^k}}\right]}_{t=0},k\ge 0,$$

(5)

and the inverse differential transforms of $x_k,y_k,k\ge 1$ are given by the following:

$$\mathbf{x}\left(t\right)=\sum _{k=0}^{\infty }{x}_k{t}^k,\mathbf{y}\left(t\right)=\sum _{k=0}^{\infty }{y}_k{t}^k.$$

(6)

From (5) and (6), we have the following:

$$\mathbf{x}\left(t\right)=\sum _{k=0}^{\infty }{\displaystyle \frac{1}{k!}}{\left[{\displaystyle \frac{d^k\mathbf{x}\left(t\right)}{d{t}^k}}\right]}_{t=0}{t}^k,\mathbf{y}\left(t\right)=\sum _{k=0}^{\infty }{\displaystyle \frac{1}{k!}}{\left[{\displaystyle \frac{d^k\mathbf{y}\left(t\right)}{d{t}^k}}\right]}_{t=0}{t}^k.$$

(7)

An approximate solution is given by the following:

$$\mathbf{x}\left(t\right)=\sum _{k=0}^K{x}_k{t}^k,\mathbf{y}\left(t\right)=\sum _{k=0}^K{y}_k{t}^k,$$

(8)

where K is the number of terms in the approximation. The accuracy of the solution should improve as we include more terms in this expansion.

If $\mathbf{x}\left(t\right)$ and $\mathbf{y}\left(t\right)$ are expanded as in (6), then the nonlinear function $G(\mathbf{x},\mathbf{y})$ can be expanded using the Adomian polynomials as follows:

$$G\left(\sum _{k=0}^{\infty }{x}_k{t}^k,\sum _{k=0}^{\infty }{y}_k{t}^k\right)=\sum _{k=0}^{\infty }{G}_k\left(x_0,x_1,\dots ,x_k;y_0,y_1,\dots ,y_k\right)t^k,$$

(9)

where $x_k$ and $y_k$, $k\ge 1$ denote the coefficients of expansion (6) and $G_k$ denotes the Adomian polynomials, as follows:

$$\begin{array}{ccc}\hfill G_0& =& G(x_0,y_0),\hfill \\ \hfill G_1& =& x_1{G}^{(1,0)}+y_1{G}^{(0,1)},\hfill \\ \hfill G_2& =& x_2{G}^{(1,0)}+y_2{G}^{(0,1)}+{\displaystyle \frac{x_1^2}{2}}{G}^{(2,0)}+{\displaystyle \frac{y_1^2}{2}}{G}^{(0,2)}+x_1{y}_1{G}^{(1,1)},\hfill \\ \hfill G_3& =& x_3{G}^{(1,0)}+y_3{G}^{(0,1)}+x_1{x}_2{G}^{(2,0)}+y_1{y}_2{G}^{(0,2)}+(x_1{y}_2+x_2{y}_1)G^{(1,1)}\hfill \\ & & +{\displaystyle \frac{x_1^3}{6}}{G}^{(3,0)}+{\displaystyle \frac{y_1^3}{6}}{G}^{(0,3)}+{\displaystyle \frac{x_1^2{y}_1}{2}}{G}^{(2,1)}+{\displaystyle \frac{x_1{y}_1^2}{2}}{G}^{(1,2)},\hfill \\ \hfill G_4& =& x_4{G}^{(1,0)}+y_4{G}^{(0,1)}+\left(x_1{x}_3+{\displaystyle \frac{x_2^2}{2}}\right)G^{(2,0)}+\left(y_1{y}_3+{\displaystyle \frac{y_2^2}{2}}\right)G^{(0,2)}\hfill \\ & & +(x_1{y}_3+x_2{y}_2+x_3{y}_1)G^{(1,1)}+{\displaystyle \frac{x_1^2{x}_2}{2}}{G}^{(3,0)}+{\displaystyle \frac{y_1^2{y}_2}{2}}{G}^{(0,3)}+\left(x_1{x}_2{y}_1+{\displaystyle \frac{x_1^2{y}_2}{2}}\right)G^{(2,1)}\hfill \\ & & +\left(x_1{y}_1{y}_2+{\displaystyle \frac{x_2{y}_1^2}{2}}\right)G^{(1,2)}+{\displaystyle \frac{x_1^4}{4!}}{G}^{(4,0)}+{\displaystyle \frac{x_1^3{y}_1}{6}}{G}^{(3,1)}+{\displaystyle \frac{x_1^2{y}_1^2}{4}}{G}^{(2,2)}\hfill \\ & & +{\displaystyle \frac{x_1{y}_1^3}{4}}{G}^{(1,3)}+{\displaystyle \frac{y_1^4}{4!}}{G}^{(0,4)},\hfill \end{array}$$

where $G^{(k,l)}:={\displaystyle \frac{{\partial }^{k+l}G\left(x_0,y_0\right)}{\partial x_0^k\partial y_0^l}},k\ge 0,l\ge 0$.

One can calculate any other Adomian polynomial from the following formula:

$$G_k:=G_k\left(x_0,x_1,\dots ,x_k;y_0,y_1,\dots ,y_k\right)={\displaystyle \frac{1}{k!}}{\displaystyle \frac{d^k}{d{\xi }^k}}{\left[G\left(\sum _{i=0}^{\infty }{\xi }^i{x}_i,\sum _{i=0}^{\infty }{\xi }^i{y}_i\right)\right]}_{\xi =0},k\ge 0.$$

(10)

One important property of the Adomian polynomials, which we will exploit to develop our algorithms, is the fact that the Adomian polynomials $G_k,k\ge 1$ are affine with respect to $x_k$ and $y_k$.

3. The New Algorithm

The proposed algorithm effectively integrates the differential transform method (DTM) with Adomian polynomials [28,29,30,31,32,33]. The primary concept is to apply the DTM directly to the system defined in class (1), where the nonlinear vector function is expanded into a power series using Adomian polynomials [34]. By leveraging the fact that Adomian polynomials $F_k$, $k\ge 1$ are affine with respect to $x_k$ and $y_k$, we derive a linear algebraic recursion system for the differential transforms of the solution. A significant benefit of this technique is that it avoids the need to convert system (1) into forms (3) or (4) before applying the DTM, greatly simplifying the algorithm. This approach is supported by two new theorems, which apply to the general first and second-order nonlinear systems given by (1).

Theorem 1.

Consider the implicit nonlinear ordinary differential system $\mathbf{F}\left({\mathbf{x}}^{\prime },\mathbf{x},t\right)=\mathbf{0}$, with $\mathbf{x}\left(0\right)={\mathbf{x}}_0\in {\mathbb{R}}^n$, $\mathbf{F}\in {\mathbb{R}}^n\times {\mathbb{R}}^n\times [0,T]$, where the Jacobian $J={\mathbf{F}}_{{\mathbf{x}}^{\prime }}({\mathbf{x}}^{\prime },\mathbf{x},t)$ is nonsingular. Assume that the function $\mathbf{F}$ is analytical with respect to its arguments and let $\mathbf{x}\left(t\right)={\sum }_{k=0}^{\infty }{x}_k{t}^k$. Let $F_k=F_k\left(x_1,\dots ,k{x}_k,(k+1)x_{k+1};x_1,\dots ,x_k\right),k\ge 1$ be the vector of k-th Adomian polynomials of the components of the vector function $\mathbf{F}$. Then, the differential transform of the solution $\mathbf{x}\left(t\right)$ is given by the linear algebraic system $kJ{x}_k=-r_{k-1}$, where $r_{k-1}=F_{k-1}(x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,x_{k-1}),k\ge 2$, with the first recursion $x_1$ being computed from $\mathbf{F}(x_1,x_0,0)=\mathbf{0}$.

Proof. 

Assume that the solution $\mathbf{x}\left(t\right)$ can be expanded as follows:

$$\mathbf{x}\left(t\right)=\sum _{k=0}^{\infty }{x}_k{t}^k,$$

(11)

where $x_k$ is the differential transform of the solution $\mathbf{x}\left(t\right)$. Then, we expand the left side of the equation $\mathbf{F}({\mathbf{x}}^{\prime }\left(t\right),\mathbf{x}\left(t\right),t)=\mathbf{0}$ in terms of the Adomian polynomials to obtain the following:

$$\mathbf{F}\left(\sum _{k=0}^{\infty }(k+1)x_{k+1}{t}^k,\sum _{k=0}^{\infty }{x}_k{t}^k,t\right)=\sum _{k=0}^{\infty }{F}_k{t}^k=0.$$

(12)

This gives the following:

$$F_k\left(x_1,\dots ,k{x}_k,(k+1)x_{k+1};x_1,\dots ,x_k\right)=\mathbf{0},k\ge 1,$$

(13)

or

$$(k+1)J{x}_{k+1}+F_k(x_1,\dots ,k{x}_k,0;x_1,\dots ,x_k),k\ge 1,$$

(14)

which leads to the following recursion system for computing the differential transform of the derivative $\mathbf{x}\left(t\right)$

$$kJ{x}_k=-F_{k-1}(x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,x_{k-1}),k\ge 2.$$

(15)

   □

Theorem 2.

Consider the implicit nonlinear ordinary differential system $\mathbf{F}({\mathbf{x}}^{″},{\mathbf{x}}^{\prime },\mathbf{x},t)=\mathbf{0}$, with $\mathbf{x}\left(0\right)={\mathbf{x}}_0,{\mathbf{x}}^{\prime }\left(0\right)={\mathbf{x}}_1\in {\mathbb{R}}^n,\mathbf{F}\in {\left({\mathbb{R}}^n\right)}^3\times [0,T]$, where the Jacobian $J={\mathbf{F}}_{{\mathbf{x}}^{″}}({\mathbf{x}}^{″},{\mathbf{x}}^{\prime },\mathbf{x},t)$ is nonsingular. Assume that the function $\mathbf{F}$ is analytical with respect to its arguments and let $\mathbf{x}\left(t\right)={\sum }_{k=0}^{\infty }{x}_k{t}^k$. Then, the differential transform of the solution $\mathbf{x}\left(t\right)$ is given by the linear algebraic system, as follows: $(k+1)(k+2)J{x}_{k+2}=-F_k\left(2x_2,\dots ,k(k+1)x_{k+1},0;2x_2,\dots ,(k+1)x_{k+1};x_1,\dots ,x_k\right),k\ge 1$, with $x_2$ being computed from $\mathbf{F}(2x_2,x_1,x_0,0)=\mathbf{0}$. Here, $F_k=F_k(6x_3,\dots ,(k+1)(k+2)x_{k+2};2x_2,\dots ,(k+1)x_{k+1};x_1,\dots ,x_k),k\ge 1$ is the vector of k-th Adomian polynomials of the components of the vector function $\mathbf{F}$.

Proof. 

Assume that the solution $\mathbf{x}\left(t\right)$ can be expanded as in (11). Then, we expand the left side of the equation $\mathbf{F}({\mathbf{x}}^{″}\left(t\right),{\mathbf{x}}^{\prime }\left(t\right),\mathbf{x}\left(t\right),t)=0$ in terms of the Adomian polynomials to obtain the following:

$$\mathbf{F}\left(\sum _{k=0}^{\infty }(k+1)(k+2)x_{k+2}{t}^k,\sum _{k=0}^{\infty }(k+1)x_{k+1}{t}^k,\sum _{k=0}^{\infty }{x}_k{t}^k,t\right)=\sum _{k=0}^{\infty }{F}_k{t}^k=\mathbf{0}.$$

(16)

This gives the following:

$$F_k\left(2x_2,\dots ,(k+1)(k+2)x_{k+2};2x_2,\dots ,(k+1)x_{k+1};x_1,\dots ,x_k\right)=0,k\ge 1,$$

(17)

or

$$(k+1)(k+2)J{x}_{k+2}+F_k(2x_2,\dots ,k(k+1)x_{k+1},0;2x_2,\dots ,(k+1)x_{k+1};x_1,\dots ,x_k),k\ge 1,$$

(18)

which leads to the following recursion system for computing the differential transform of the solution $\mathbf{x}\left(t\right)$

$$(k-1)kJ{x}_k=-F_{k-2}\left(2x_2,\dots ,(k-2)(k-1)x_{k-1},0;2x_2,\dots ,(k-1)x_{k-1};x_1,\dots ,x_{k-2}\right),k\ge 3.$$

(19)

   □

Finally, we have Algorithms 1–3:   

Algorithm 1: DTM algorithm for solution of $\mathbf{F}({\mathbf{x}}^{\prime },\mathbf{x},t)=0,\mathbf{x}\left(0\right)={\mathbf{x}}_0$.

Input  $K,n,{\mathbf{x}}_0\in {\mathbb{R}}^n,\mathbf{F}\in {\mathbb{R}}^n$ Output  approximate solution: ${\displaystyle \mathbf{x}\left(t\right)=\sum _{k=0}^K{x}_k{t}^k}$ For $2\le k\le K$ do compute $r_{k-1}:=F_{k-1}(x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,x_{k-1})$ end do

Initialization $x_0:={\mathbf{x}}_0$ solve for $x_1$ the algebraic system: $\mathbf{F}(x_1,x_0,0)=\mathbf{0}$ compute $J={\mathbf{F}}_{{\mathbf{x}}^{\prime }}(x_1,x_0,0)$ For $2\le k\le K$ do solve for $x_k$ the linear algebraic system: $kJ{x}_k=-r_{k-1}$ end do

Algorithm 2: DTM algorithm for solution of $\mathbf{F}({\mathbf{x}}^{″},{\mathbf{x}}^{\prime },\mathbf{x},t)=\mathbf{0},\mathbf{x}\left(0\right)={\mathbf{x}}_0,{\mathbf{x}}^{\prime }\left(0\right)={\mathbf{x}}_1$.

Input  $K,n,{\mathbf{x}}_0,{\mathbf{x}}_1\in {\mathbb{R}}^n,\mathbf{F}\in {\mathbb{R}}^n$ Output  approximate solution: ${\displaystyle \mathbf{x}\left(t\right)=\sum _{k=0}^K{x}_k{t}^k}$ For $3\le k\le K$ do compute $r_{k-2}:=F_{k-2}(2x_2,\dots ,(k-2)(k-1)x_{k-1},0;2x_2,\dots ,(k-1)x_{k-1};x_1,\dots ,x_{k-2})$ end do

Initialization $x_0:={\mathbf{x}}_0,x_1:={\mathbf{x}}_1$ solve for $x_2$ the algebraic system: $\mathbf{F}(2x_2,x_1,x_0,0)=\mathbf{0}$ compute $J={\mathbf{F}}_{{\mathbf{x}}^{″}}(2x_2,x_1,x_0,0)$ For $3\le k\le K$ do solve for $x_k$ the linear algebraic system: $k(k-1)J{x}_k=-r_{k-2}$ end do

We can extend the interval of convergence of the solution by using a multi-stage DTM (MsDTM). The procedure of the MsDTM for $\mathbf{F}\left({\mathbf{x}}^{\prime },\mathbf{x},t\right)=\mathbf{0}$, with $\mathbf{x}\left(0\right)={\mathbf{x}}_0$ is described in the following algorithm:

Algorithm 3: MsDTM algorithm.

Choose the number N of sub-divisions and the order of approximation K. Sub-divide the interval $\left[0,T\right]$ into N equal sub-intervals $\left[t_{i-1},t_i\right]$, of size $h=T/N$, where $t_i=ih$, $i=1,\dots ,N$. For $i=1$, apply the DTM to $\mathbf{F}\left({\mathbf{x}}^{\prime },\mathbf{x},t\right)=\mathbf{0}$ on $\left[0,t_i\right]$ using ${\mathbf{x}}_0^i=\mathbf{x}\left(0\right)$ to obtain an approximate solution ${\mathbf{x}}^i\left(t\right)={\sum }_{k=0}^K{x}_k^i{t}^k$. From $i=2$ to $i=N$, apply the DTM to $\mathbf{F}\left({\mathbf{x}}^{\prime },\mathbf{x},t\right)=\mathbf{0}$ on $\left[t_{i-1},t_i\right]$ using ${\mathbf{x}}_0^i={\mathbf{x}}^{i-1}\left(t_{i-1}\right)$ to obtain an approximate solution ${\mathbf{x}}^i\left(t\right)={\sum }_{k=0}^K{x}_k^i{(t-t_{i-1})}^k$. Repeat step 4. An approximate solution to the initial-value problem $\mathbf{F}\left({\mathbf{x}}^{\prime },\mathbf{x},t\right)=\mathbf{0}$, $\mathbf{x}\left(0\right)={\mathbf{x}}_0$ on $[0,T]$ is given by $\mathbf{x}\left(t\right)=\left\{\begin{array}{cc}{\mathbf{x}}^1\left(t\right),\hfill & 0\le t\le t_1\hfill \\ {\mathbf{x}}^2\left(t\right),\hfill & t_1<t\le t_2\hfill \\ ⋮\hfill & \\ {\mathbf{x}}^N\left(t\right),\hfill & t_{N-1}<t\le T\hfill \end{array}\right.$

A similar multi-stage algorithm can be given for second-order implicit systems $\mathbf{F}\left({\mathbf{x}}^{″},{\mathbf{x}}^{\prime },\mathbf{x},t\right)=\mathbf{0}$.

4. Numerical Results

In this section, we showcase the effectiveness and accuracy of our proposed technique by solving several first and second-order implicit ODEs and systems of ODEs with different nonlinearities. Notably, these examples cannot be solved using Maple 2024, either numerically or through power series expansions. Additionally, conventional methods like Runge–Kutta do not apply to these types of ODEs and systems. Utilizing Algorithms 1 and 2, our technique successfully computes the solutions in power series form. The implementation of our algorithms was carried out in Maple 2024. Our method is capable of providing the exact power series expansion of the solution, provided that all computations are performed with precision. To solve these examples over larger intervals, we used the multi-stage versions of Algorithms 1 and 2. The numerical results show that the MsDTM provides accurate solutions to these implicit systems of ODEs.

Example 1.

In this example, we consider the following fully implicit nonlinear system of first-order ordinary differential equations, as follows:

$$\begin{array}{ccc}\hfill {\mathbf{x}}_1^{\prime }-{\displaystyle \frac{{\mathbf{x}}_1}{{{\mathbf{x}}^{\prime }}_2^2+{\mathbf{x}}_2^2}}+{\mathbf{x}}_1-{\mathbf{x}}_2& =& 0,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\mathbf{x}}_2^{\prime }-{\displaystyle \frac{{\mathbf{x}}_2}{{{\mathbf{x}}^{\prime }}_1^2+{\mathbf{x}}_1^2}}+{\mathbf{x}}_1+{\mathbf{x}}_2& =& 0,\hfill \end{array}$$

(21)

where $t>0$. The sought solution is as follows: $\mathbf{x}={\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)}^T$. The differential system (20) and (21) has the form (1), with $n=2$, and $m=1$. This system is subject to the following initial conditions:

$$\mathbf{x}\left(0\right)=\left(\begin{array}{c}0\\ 1\end{array}\right).$$

(22)

We should emphasize here that Maple 2024 software cannot solve this system. However, our technique has successfully obtained the exact solution of this implicit system in power series. By applying Algorithm 1 to (20) and (21) with initial condition (22), we can recursively determine the DTM expansion coefficients $x_k$ for $k\ge 2$. According to the initial condition (22), we have the following:

$$x_0=\left(\begin{array}{c}0\\ 1\end{array}\right).$$

(23)

Then, we determine the differential transform $x_1={\left(x_{11},x_{12}\right)}^T$ from the algebraic system as follows:

$$\mathbf{F}(x_1,x_0,0)=\mathbf{0}.$$

(24)

The above system has a unique solution $x_1$ and simplifies to the following:

$$\begin{array}{ccc}\hfill x_{11}-{\displaystyle \frac{x_{01}}{x_{12}^2+x_{02}^2}}+x_{01}-x_{02}& =& 0,\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill x_{12}-{\displaystyle \frac{x_{02}}{x_{11}^2+x_{01}^2}}+x_{01}+x_{02}& =& 0.\hfill \end{array}$$

(26)

Solving this algebraic system, we obtain the following:

$$x_1=\left(\begin{array}{c}{x}_{11}\\ x_{12}\end{array}\right)=\left(\begin{array}{c}1\\ 0\end{array}\right).$$

(27)

Then, we compute the Jacobian $J:={\mathbf{F}}_{{\mathbf{x}}^{\prime }}$ of the function $\mathbf{F}$ with respect to ${\mathbf{x}}^{\prime }$ at the point $(x_1,x_0,0)$. We obtain the following:

$$J=\left(\begin{array}{cc}1& 0\\ 2& 1\end{array}\right).$$

(28)

Using Theorem 1, we have the following recursion for the differential transform $x_k$

$$x_k=-{\left(kJ\right)}^{-1}{r}_{k-1},k\ge 2,$$

(29)

where

$$r_{k-1}:=F_{k-1}(x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,\dots ,x_{k-1}).$$

(30)

Using Algorithm 1, we compute the inverse of the matrix and the right-hand sides $r_{k-1}$ for $k\ge 2$ and deduce the following:

For $k=2$, we have the following:

$$x_2=-{\left(2J\right)}^{-1}{r}_1=-\left(\begin{array}{cc}{\displaystyle \frac{1}{2}}& 0\\ -1& {\displaystyle \frac{1}{2}}\end{array}\right)\left(\begin{array}{c}0\\ 1\end{array}\right)=\left(\begin{array}{c}0\\ -{\displaystyle \frac{1}{2}}\end{array}\right).$$

(31)

For $k=3$, we have the following:

$$x_3=-{\left(3J\right)}^{-1}{r}_2=-\left(\begin{array}{cc}{\displaystyle \frac{1}{3}}& 0\\ -{\displaystyle \frac{2}{3}}& {\displaystyle \frac{1}{3}}\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{1}{2}}\\ 1\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{6}}\\ 0\end{array}\right).$$

(32)

For $k=4$, we have the following:

$$x_4=-{\left(4J\right)}^{-1}{r}_3=-\left(\begin{array}{cc}{\displaystyle \frac{1}{4}}& 0\\ -{\displaystyle \frac{1}{2}}& {\displaystyle \frac{1}{4}}\end{array}\right)\left(\begin{array}{c}0\\ -{\displaystyle \frac{1}{6}}\end{array}\right)=\left(\begin{array}{c}0\\ {\displaystyle \frac{1}{24}}\end{array}\right).$$

(33)

For $k=5$, we have the following:

$$x_5=-{\left(5J\right)}^{-1}{r}_4=-\left(\begin{array}{cc}{\displaystyle \frac{1}{5}}& 0\\ -{\displaystyle \frac{2}{5}}& {\displaystyle \frac{1}{5}}\end{array}\right)\left(\begin{array}{c}-{\displaystyle \frac{1}{24}}\\ -{\displaystyle \frac{1}{12}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{120}}\\ 0\end{array}\right).$$

(34)

For $k=6$, we have the following:

$$x_6=-{\left(6J\right)}^{-1}{r}_5=-\left(\begin{array}{cc}{\displaystyle \frac{1}{6}}& 0\\ -{\displaystyle \frac{1}{3}}& {\displaystyle \frac{1}{6}}\end{array}\right)\left(\begin{array}{c}0\\ {\displaystyle \frac{1}{120}}\end{array}\right)=\left(\begin{array}{c}0\\ -{\displaystyle \frac{1}{720}}\end{array}\right).$$

(35)

Finally, we construct the following approximate solution with $K=7$:

$$\mathbf{x}\left(t\right)=\left(\begin{array}{c}0\\ 1\end{array}\right)+t\left(\begin{array}{c}1\\ 0\end{array}\right)+t^2\left(\begin{array}{c}0\\ -{\displaystyle \frac{1}{2}}\end{array}\right)+t^3\left(\begin{array}{c}-{\displaystyle \frac{1}{6}}\\ 0\end{array}\right)+t^4\left(\begin{array}{c}0\\ {\displaystyle \frac{1}{24}}\end{array}\right)+t^5\left(\begin{array}{c}{\displaystyle \frac{1}{120}}\\ 0\end{array}\right)+t^6\left(\begin{array}{c}0\\ -{\displaystyle \frac{1}{720}}\end{array}\right).$$

(36)

One can easily check that the above approximation forms the first terms of the Maclaurin power series of the exact solution $\mathbf{x}\left(t\right)={(\sin \left(t\right),\cos \left(t\right))}^T$. It is worth noting here that we obtain the exact values of the successive DTM expansion coefficients. It is known that by increasing the number of terms in the series solution, the error decreases, and the convergence interval of the truncated series increases. To compute accurate numerical solutions over large internals and extend the interval of convergence of the solution, we apply the multi-stage DTM (MsDTM) algorithm given in Section 3. To solve this example with the MsDTM, we use an order of approximation $K=7$ and subdivide the solution interval $[0,T]=[0,15]$ into $N=400$ sub-intervals. Starting the recursion with $x_0=x\left(0\right)={\left(0,1\right)}^T$, we obtained the numerical results shown in Figure 1 and Figure 2. Each of these figures shows the numerical solution (solid line), the exact solution (dashed line), and the corresponding absolute error. From the graph (a), we can see that the numerical solution approximates well the exact solution. In each graph (b), we can see the corresponding absolute error. These graphs show that the errors of approximation using the MsDTM are less than $1\times {10}^{-11}$. These numerical results show that equations in this example are satisfied with high accuracy and that the MsDTM computes accurate numerical solutions over large intervals. To increase the accuracy in the approximate MsDTM solution, one can increase the approximation order K or the number N of sub-divisions of the solution interval $[0,T]$. To maintain a given accuracy while increasing the size of the solution interval, one may need to increase the number of terms K and or increase the number N of interval subdivisions.

Example 2.

In this example, we consider the following fully implicit nonlinear system of first-order ordinary differential equations:

$$\begin{array}{ccc}\hfill 10{\mathbf{x}}_1^{\prime }-\cos \left(\pi {\mathbf{x}}_2^{\prime }\right)+{\displaystyle \frac{5}{3}}{\mathbf{x}}_1+10{\mathbf{x}}_2+\cos \left(\pi {\mathbf{x}}_1-{\displaystyle \frac{\pi }{6}}{\mathbf{x}}_2\right)& =& 0,\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill 10{\mathbf{x}}_2^{\prime }+\sin \left(\pi {\mathbf{x}}_1^{\prime }\right)-10{\mathbf{x}}_1+{\displaystyle \frac{5}{3}}{\mathbf{x}}_2+\sin \left({\displaystyle \frac{\pi }{6}}{\mathbf{x}}_1+\pi {\mathbf{x}}_2\right)& =& 0,\hfill \end{array}$$

(38)

where $t>0$. The sought solution is as follows: $\mathbf{x}={\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)}^T$. The differential system (37) and (38) has the form (1), with $n=2$, and $m=1$. This system is subject to the following initial conditions:

$$\mathbf{x}\left(0\right)=\left(\begin{array}{c}1\\ 0\end{array}\right).$$

(39)

We should emphasize here that Maple 2024 software cannot solve this system. However, our technique has successfully obtained the exact solution of this implicit system in power series. By applying Algorithm 1 to (37) and (38) with initial condition (39), we can recursively determine the DTM expansion coefficients $x_k$ for $k\ge 2$. According to the initial condition (39), we have the following:

$$x_0=\left(\begin{array}{c}1\\ 0\end{array}\right).$$

(40)

Then, we determine the differential transform $x_1={\left(x_{11},x_{12}\right)}^T$ from the algebraic system, as follows:

$$\mathbf{F}(x_1,x_0,0)=\mathbf{0}.$$

(41)

The above system has a unique solution $x_1$ and simplifies to the following:

$$\begin{array}{ccc}\hfill 10x_{11}-\cos \left(\pi x_{12}\right)+{\displaystyle \frac{5}{3}}+\cos \left(\pi \right)& =& 0,\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill 10x_{12}+\sin \left(\pi x_{11}\right)-10+\sin \left({\displaystyle \frac{\pi }{6}}\right)& =& 0.\hfill \end{array}$$

(43)

Solving this algebraic system, we obtain the following:

$$x_1=\left(\begin{array}{c}{x}_{11}\\ x_{12}\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{6}}\\ 1\end{array}\right).$$

(44)

Then, we compute the Jacobian $J:={\mathbf{F}}_{{\mathbf{x}}^{\prime }}$ of the function $\mathbf{F}$ with respect to ${\mathbf{x}}^{\prime }$ at the point $(x_1,x_0,0)$. We obtain the following:

$$J=\left(\begin{array}{cc}{\displaystyle \frac{1}{2}}\pi \sqrt{3}& 10\\ 10& 0\end{array}\right).$$

(45)

Using Theorem 1, we have the following recursion for the differential transform $x_k$:

$$x_k=-{\left(kJ\right)}^{-1}{r}_{k-1},k\ge 2,$$

(46)

where

$$r_{k-1}:=F_{k-1}(x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,\dots ,x_{k-1}).$$

(47)

Using Algorithm 1, we compute the inverse of the matrix and the right-hand sides $r_{k-1}$ for $k\ge 2$ and deduce the following:

For $k=2$, we have the following:

$$x_2=-{\left(2J\right)}^{-1}{r}_1=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{20}}\\ {\displaystyle \frac{1}{20}}& -{\displaystyle \frac{1}{400}}\pi \sqrt{3}\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{10}{3}}+{\displaystyle \frac{35}{72}}\pi \sqrt{3}\\ {\displaystyle \frac{175}{18}}\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{35}{72}}\\ -{\displaystyle \frac{1}{6}}\end{array}\right).$$

(48)

For $k=3$, we have the following:

$$x_3=-{\left(3J\right)}^{-1}{r}_2=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{30}}\\ {\displaystyle \frac{1}{30}}& -{\displaystyle \frac{1}{600}}\pi \sqrt{3}\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{55}{12}}-{\displaystyle \frac{107}{864}}\pi \sqrt{3}\\ -{\displaystyle \frac{535}{216}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{107}{1296}}\\ -{\displaystyle \frac{11}{72}}\end{array}\right).$$

(49)

For $k=4$, we have the following:

$$x_4=-{\left(4J\right)}^{-1}{r}_3=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{40}}\\ {\displaystyle \frac{1}{40}}& -{\displaystyle \frac{1}{800}}\pi \sqrt{3}\end{array}\right)\left(\begin{array}{c}-{\displaystyle \frac{175}{162}}-{\displaystyle \frac{1081}{15552}}\pi \sqrt{3}\\ -{\displaystyle \frac{5405}{3888}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1081}{31104}}\\ {\displaystyle \frac{35}{1296}}\end{array}\right).$$

(50)

For $k=5$, we have the following:

$$x_5=-{\left(5J\right)}^{-1}{r}_4=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{50}}\\ {\displaystyle \frac{1}{50}}& -{\displaystyle \frac{1}{1000}}\pi \sqrt{3}\end{array}\right)\left(\begin{array}{c}-{\displaystyle \frac{4705}{15552}}+{\displaystyle \frac{6121}{373248}}\pi \sqrt{3}\\ {\displaystyle \frac{30605}{93312}}\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{6121}{933120}}\\ {\displaystyle \frac{941}{155520}}\end{array}\right).$$

(51)

For $k=6$, we have the following:

$$x_6=-{\left(6J\right)}^{-1}{r}_5=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{60}}\\ {\displaystyle \frac{1}{60}}& -{\displaystyle \frac{1}{1200}}\pi \sqrt{3}\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{1177}{15552}}+{\displaystyle \frac{5551}{2239488}}\pi \sqrt{3}\\ {\displaystyle \frac{27755}{559872}}\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{5551}{6718464}}\\ -{\displaystyle \frac{1177}{933120}}\end{array}\right).$$

(52)

Finally, we construct the following approximate solution with $K=7$:

$$\mathbf{x}\left(t\right)=\left(\begin{array}{c}1\\ 0\end{array}\right)+t\left(\begin{array}{c}-{\displaystyle \frac{1}{6}}\\ 1\end{array}\right)+t^2\left(\begin{array}{c}-{\displaystyle \frac{35}{72}}\\ -{\displaystyle \frac{1}{6}}\end{array}\right)+t^3\left(\begin{array}{c}{\displaystyle \frac{107}{1296}}\\ -{\displaystyle \frac{11}{72}}\end{array}\right)+t^4\left(\begin{array}{c}{\displaystyle \frac{1081}{31104}}\\ {\displaystyle \frac{35}{1296}}\end{array}\right)+t^5\left(\begin{array}{c}-{\displaystyle \frac{6121}{933120}}\\ {\displaystyle \frac{941}{155520}}\end{array}\right)+t^6\left(\begin{array}{c}-{\displaystyle \frac{5551}{6718464}}\\ -{\displaystyle \frac{1177}{933120}}\end{array}\right).$$

(53)

One can easily check that the above approximation forms the first terms of the Maclaurin power series of the exact solution $\mathbf{x}\left(t\right)={(e^{-t/6}\cos \left(t\right),e^{-t/6}\sin \left(t\right))}^T$. It is worth noting here that we obtain the exact values of the successive DTM expansion coefficients. It is known that by increasing the number of terms in the series solution, the error decreases, and the convergence interval of the truncated series increases (Figure 3 and Figure 4).

Example 3.

In this example, we consider the following fully implicit nonlinear system of first-order ordinary differential equations:

$$\begin{array}{ccc}\hfill \arctan \left(t{\mathbf{x}}_1^{\prime }\right)-{\mathbf{x}}_2^{\prime }-{\displaystyle \frac{1}{1+t^2}}-\arctan \left({\displaystyle \frac{t}{1+t^2}}\right)& =& 0,\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill {\mathbf{x}}_1^{\prime }+\arctan \left(t{\mathbf{x}}_2^{\prime }+{\displaystyle \frac{t}{1+t^2}}\right)-{\displaystyle \frac{1}{1+t^2}}& =& 0,\hfill \end{array}$$

(55)

where $t>0$. The sought solution is $\mathbf{x}={\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)}^T$. The differential system (54) and (55) has the form (1), with $n=2$, and $m=1$. This system is subject to the following initial conditions:

$$\mathbf{x}\left(0\right)=\left(\begin{array}{c}1\\ 1\end{array}\right).$$

(56)

We should emphasize here that Maple 2024 software cannot solve this system. However, our technique has successfully obtained the exact solution of this implicit system in power series. By applying Algorithm 1 to (54) and (55) with initial condition (56), we can recursively determine the DTM expansion coefficients $x_k$ for $k\ge 2$. According to the initial condition (56), we have the following:

$$x_0=\left(\begin{array}{c}1\\ 1\end{array}\right).$$

(57)

Then, we determine the differential transform $x_1={\left(x_{11},x_{12}\right)}^T$ from the algebraic system as follows:

$$\mathbf{F}(x_1,x_0,0)=\mathbf{0}.$$

(58)

The above system has a unique solution $x_1$ and simplifies to the following:

$$\begin{array}{ccc}\hfill x_{12}+1& =& 0,\hfill \end{array}$$

(59)

$$\begin{array}{ccc}\hfill x_{11}-1& =& 0.\hfill \end{array}$$

(60)

Solving this algebraic system, we obtain the following:

$$x_1=\left(\begin{array}{c}{x}_{11}\\ x_{12}\end{array}\right)=\left(\begin{array}{c}1\\ -1\end{array}\right).$$

(61)

Then, we compute the Jacobian $J:={\mathbf{F}}_{{\mathbf{x}}^{\prime }}$ of the function $\mathbf{F}$ with respect to ${\mathbf{x}}^{\prime }$ at the point $(x_1,x_0,0)$. We obtain the following:

$$J=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right).$$

(62)

Using Theorem 1, we have the following recursion for the differential transform: $x_k$

$$x_k=-{\left(kJ\right)}^{-1}{r}_{k-1},k\ge 2,$$

(63)

where

$$r_{k-1}:=F_{k-1}(x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,(k-1)x_{k-1},0;x_1,\dots ,\dots ,x_{k-1}).$$

(64)

Using Algorithm 1, we compute the inverse of the matrix and the right-hand sides $r_{k-1}$ for $k\ge 2$ and deduce the following:

For $k=2$, we have the following:

$$x_2=-{\left(2J\right)}^{-1}{r}_1=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{2}}\\ -{\displaystyle \frac{1}{2}}& 0\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(65)

For $k=3$, we have the following:

$$x_3=-{\left(3J\right)}^{-1}{r}_2=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{3}}\\ -{\displaystyle \frac{1}{3}}& 0\end{array}\right)\left(\begin{array}{c}1\\ 1\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{3}}\\ {\displaystyle \frac{1}{3}}\end{array}\right).$$

(66)

For $k=4$, we have the following:

$$x_4=-{\left(4J\right)}^{-1}{r}_3=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{4}}\\ -{\displaystyle \frac{1}{4}}& 0\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(67)

For $k=5$, we have the following:

$$x_5=-{\left(5J\right)}^{-1}{r}_4=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{5}}\\ -{\displaystyle \frac{1}{5}}\end{array}\right)\left(\begin{array}{c}-1\\ -1\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{5}}\\ -{\displaystyle \frac{1}{5}}\end{array}\right).$$

(68)

For $k=6$, we have the following:

$$x_6=-{\left(6J\right)}^{-1}{r}_5=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{6}}\\ -{\displaystyle \frac{1}{6}}& 0\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(69)

For $k=7$, we have the following:

$$x_7=-{\left(7J\right)}^{-1}{r}_6=-\left(\begin{array}{cc}0& {\displaystyle \frac{1}{7}}\\ -{\displaystyle \frac{1}{7}}& 0\end{array}\right)\left(\begin{array}{c}1\\ 1\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{7}}\\ {\displaystyle \frac{1}{7}}\end{array}\right).$$

(70)

Finally, we construct the following approximate solution with $K=7$:

$$\mathbf{x}\left(t\right)=\left(\begin{array}{c}1\\ 1\end{array}\right)+t\left(\begin{array}{c}1\\ -1\end{array}\right)+t^3\left(\begin{array}{c}-{\displaystyle \frac{1}{3}}\\ {\displaystyle \frac{1}{3}}\end{array}\right)+t^5\left(\begin{array}{c}{\displaystyle \frac{1}{5}}\\ -{\displaystyle \frac{1}{5}}\end{array}\right)+t^7\left(\begin{array}{c}-{\displaystyle \frac{1}{7}}\\ {\displaystyle \frac{1}{7}}\end{array}\right).$$

(71)

One can easily check that the above approximation forms the first terms of the Maclaurin power series of the exact solution $\mathbf{x}\left(t\right)={(1+\arctan \left(t\right),1-\arctan \left(t\right))}^T$. It is worth noting here that we obtain the exact values of the successive DTM expansion coefficients. It is known that by increasing the number of terms in the series solution, the error decreases, and the convergence interval of the truncated series increases (Figure 5 and Figure 6).

Example 4.

In this example, we consider the following second-order implicit nonlinear equation:

$$2{\mathbf{x}}^{″}\left(t\right)+\sin ({\mathbf{x}}^{″}\left(t\right)-t^3)+4{\mathbf{x}}^{\prime }\left(t\right)-t^4-2t^3-4=0,t\ge 0,$$

(72)

subject to the initial conditions

$$\mathbf{x}\left(0\right)=0,{\mathbf{x}}^{\prime }\left(0\right)=1.$$

(73)

We should emphasize here that Maple 2024 software could not solve this equation numerically or in power series in its present form. However, there is no trouble for our technique to obtain the solution of this implicit equation in power series. By applying Algorithm 2 to (72) with initial condition (73), we can recursively determine the DTM expansion coefficients $x_k$ for $k\ge 2$. The differential Equation (72) has the following form:

$$\mathbf{F}({\mathbf{x}}^{″},{\mathbf{x}}^{\prime },\mathbf{x},t)=0,$$

(74)

where $\mathbf{F}({\mathbf{x}}^{″},{\mathbf{x}}^{\prime },\mathbf{x},t)=2{\mathbf{x}}^{″}+\sin ({\mathbf{x}}^{″}-t^3)+4{\mathbf{x}}^{\prime }-t^4-2t^3-4$. According to the initial conditions (73), we have the following:

$$x_0=\mathbf{x}\left(0\right)=0,x_1={\mathbf{x}}^{\prime }\left(0\right)=1.$$

(75)

Then, we determine the second DTM expansion coefficient $x_2$ from the algebraic equation, as follows:

$$\mathbf{F}(2x_2,x_1,x_0,0)=0.$$

(76)

This equation has a unique solution $x_2$ and simplifies to the following:

$$4x_2+\sin \left(2x_2\right)=0.$$

(77)

Solving this equation, we obtain the following:

$$x_2=0.$$

(78)

Then, we compute the Jacobian J of the function $\mathbf{F}$ with respect to ${\mathbf{x}}^{″}$ at the point $(2x_2,x_1,x_0,0)$. For the present example, J is just a scalar, as follows:

$$J=2+\cos \left(2x_2\right)=3.$$

(79)

Using Theorem 2, we have the following recursion for the DTM expansion coefficients $x_k$:

$$x_k=-\left({\displaystyle \frac{1}{3k(k-1)}}\right)r_{k-2},k\ge 3,$$

(80)

where

$$r_{k-2}:=F_{k-2}(2x_2,\dots ,(k-2)(k-1)x_{k-1},0;2x_2,\dots ,(k-1)x_{k-1};x_1,\dots ,x_{k-2}).$$

(81)

Using Algorithm 2, we compute the right-hand sides $r_{k-2}$ for $k\ge 3$ and obtain the following:

$$r_1=8x_2=0.$$

(82)

Using (80), we deduce the following:

$$x_3=0.$$

(83)

Then, we have the following:

$$r_2=12x_3-18x_3^2\sin \left(2x_2\right)=0.$$

(84)

Using (80), we deduce the following:

$$x_4=0.$$

(85)

Then, we have the following:

$$r_3=-2+16x_4-\cos \left(2x_2\right)-72x_4{x}_3\sin \left(2x_2\right)-36x_3^3\cos \left(2x_2\right)=-3.$$

(86)

Using (80), we deduce the following:

$$x_5={\displaystyle \frac{1}{20}}.$$

(87)

In a similar way, we can show that $r_k=0$ for $k\ge 4$, hence, $x_k=0$ for $k\ge 6$. Finally, we construct the following approximate solution based on five terms:

$$\mathbf{x}\left(t\right)=\sum _{k=0}^5{x}_k{t}^k=t+{\displaystyle \frac{1}{20}}{t}^5.$$

(88)

One can easily check that the exact solution is $\mathbf{x}\left(t\right)=t+{\displaystyle \frac{1}{20}}{t}^5$. It is worth noting here that we obtain the exact values of the successive DTM expansion coefficients.

Example 5.

In this example, we consider the following fully implicit nonlinear system of second-order ordinary differential equations:

$$\begin{array}{ccc}\hfill 3{\mathbf{x}}_{\mathbf{1}}^{″}+{\mathbf{x}}_{\mathbf{2}}^{″}+\sin ({\mathbf{x}}_{\mathbf{1}}^{″}+2{\mathbf{x}}_1-2{\mathbf{x}}_2)-{\mathbf{x}}_{\mathbf{2}}^{\prime 2}+4({\mathbf{x}}_1-{\mathbf{x}}_2+{\mathbf{x}}_1{\mathbf{x}}_2)& =& 0,\hfill \end{array}$$

(89)

$$\begin{array}{ccc}\hfill 4{\mathbf{x}}_{\mathbf{2}}^{″}+\sin ({\mathbf{x}}_{\mathbf{1}}^{″}-2\cos \left(2t\right))-{\mathbf{x}}_{\mathbf{1}}^{\prime 2}+4(2{\mathbf{x}}_{\mathbf{2}}-2{\mathbf{x}}_1+{\mathbf{x}}_1{\mathbf{x}}_2)& =& 0,\hfill \end{array}$$

(90)

where $t>0$. The sought solution is as follows: $\mathbf{x}={\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)}^T$. The differential system (89) and (90) has the form (1), with $n=2$, and $m=2$. This system is subject to the following initial conditions:

$$\mathbf{x}\left(0\right)=\left(\begin{array}{c}0\\ 1\end{array}\right),{\mathbf{x}}^{\prime }\left(0\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(91)

We should emphasize here that Maple 2024 software cannot solve this system. However, our technique has successfully obtained the exact solution of this implicit system in power series. By applying Algorithm 2 to (89) and (90) with initial condition (91), we can recursively determine the DTM expansion coefficients $x_k$ for $k\ge 3$. According to the initial condition (91), we have the following:

$$x_0=\left(\begin{array}{c}0\\ 1\end{array}\right),x_1=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(92)

Then, we determine the differential transform $x_2={\left(x_{21},x_{22}\right)}^T$ from the algebraic system, as follows:

$$\mathbf{F}(2x_2,x_1,x_0,0)=\mathbf{0}.$$

(93)

The above system has a unique solution $x_2$ and simplifies to the following:

$$\begin{array}{ccc}\hfill 6x_{21}+2x_{22}+\sin (2x_{21}-2)-4& =& 0,\hfill \end{array}$$

(94)

$$\begin{array}{ccc}\hfill 8x_{22}+\sin (2x_{21}-2)+8& =& 0.\hfill \end{array}$$

(95)

Solving this algebraic system, we obtain the following:

$$x_2=\left(\begin{array}{c}{x}_{21}\\ x_{22}\end{array}\right)=\left(\begin{array}{c}1\\ -1\end{array}\right).$$

(96)

Then, we compute the Jacobian $J:={\mathbf{F}}_{{\mathbf{x}}^{″}}$ of the function $\mathbf{F}$ with respect to ${\mathbf{x}}^{″}$ at the point $(x_2,x_1,x_0,0)$. We obtain the following:

$$J=\left(\begin{array}{cc}4& 1\\ 1& 4\end{array}\right).$$

(97)

Using Theorem 2, we have the following recursion for the differential transform: $x_k$

$$x_k=-{\left(k(k-1)J\right)}^{-1}{r}_{k-2},k\ge 3,$$

(98)

where

$$r_{k-2}:=F_{k-2}(2x_2,\dots ,(k-2)(k-1)x_{k-1},0;2x_2,\dots ,(k-1)x_{k-1};x_1,\dots ,\dots ,x_{k-2}).$$

(99)

Using Algorithm 2, we compute the inverse of the matrix and the right-hand sides $r_{k-2}$ for $k\ge 3$ and deduce the following:

For $k=3$, we have the following:

$$x_3=-{\left(6J\right)}^{-1}{r}_1=\left(\begin{array}{cc}-{\displaystyle \frac{2}{45}}& {\displaystyle \frac{1}{90}}\\ {\displaystyle \frac{1}{90}}& -{\displaystyle \frac{2}{45}}\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(100)

For $k=4$, we have the following:

$$x_4=-{\left(12J\right)}^{-1}{r}_2=\left(\begin{array}{cc}-{\displaystyle \frac{1}{45}}& {\displaystyle \frac{1}{180}}\\ {\displaystyle \frac{1}{180}}& -{\displaystyle \frac{1}{45}}\end{array}\right)\left(\begin{array}{c}12\\ -12\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{3}}\\ {\displaystyle \frac{1}{3}}\end{array}\right).$$

(101)

For $k=5$, we have the following:

$$x_5=-{\left(20J\right)}^{-1}{r}_3=\left(\begin{array}{cc}-{\displaystyle \frac{1}{75}}& {\displaystyle \frac{1}{300}}\\ {\displaystyle \frac{1}{300}}& -{\displaystyle \frac{1}{75}}\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(102)

For $k=6$, we have the following:

$$x_6=-{\left(30J\right)}^{-1}{r}_4=\left(\begin{array}{cc}-{\displaystyle \frac{2}{225}}& {\displaystyle \frac{1}{450}}\\ {\displaystyle \frac{1}{450}}& -{\displaystyle \frac{2}{225}}\end{array}\right)\left(\begin{array}{c}-4\\ 4\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{2}{45}}\\ -{\displaystyle \frac{2}{45}}\end{array}\right).$$

(103)

For $k=7$, we have the following:

$$x_7=-{\left(42J\right)}^{-1}{r}_5=\left(\begin{array}{cc}-{\displaystyle \frac{2}{315}}& {\displaystyle \frac{1}{630}}\\ {\displaystyle \frac{1}{630}}& -{\displaystyle \frac{2}{315}}\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(104)

For $k=8$, we have the following:

$$x_8=-{\left(56J\right)}^{-1}{r}_6=\left(\begin{array}{cc}-{\displaystyle \frac{1}{210}}& {\displaystyle \frac{1}{840}}\\ {\displaystyle \frac{1}{840}}& -{\displaystyle \frac{1}{210}}\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{8}{15}}\\ -{\displaystyle \frac{8}{15}}\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{315}}\\ {\displaystyle \frac{1}{315}}\end{array}\right).$$

(105)

For $k=9$, we have the following:

$$x_9=-{\left(72J\right)}^{-1}{r}_7=\left(\begin{array}{cc}-{\displaystyle \frac{1}{270}}& {\displaystyle \frac{1}{1080}}\\ {\displaystyle \frac{1}{1080}}& -{\displaystyle \frac{1}{270}}\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(106)

For $k=10$, we have the following:

$$x_{10}=-{\left(90J\right)}^{-1}{r}_8=\left(\begin{array}{cc}-{\displaystyle \frac{2}{675}}& {\displaystyle \frac{1}{1350}}\\ {\displaystyle \frac{1}{1350}}& -{\displaystyle \frac{2}{675}}\end{array}\right)\left(\begin{array}{c}-{\displaystyle \frac{4}{105}}\\ {\displaystyle \frac{4}{105}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{2}{14175}}\\ -{\displaystyle \frac{2}{14175}}\end{array}\right).$$

(107)

Finally, we construct the following approximate solution with $K=10$:

$$\mathbf{x}\left(t\right)=\left(\begin{array}{c}0\\ 1\end{array}\right)+t^2\left(\begin{array}{c}1\\ -1\end{array}\right)+t^4\left(\begin{array}{c}-{\displaystyle \frac{1}{3}}\\ {\displaystyle \frac{1}{3}}\end{array}\right)+t^6\left(\begin{array}{c}{\displaystyle \frac{2}{45}}\\ -{\displaystyle \frac{2}{45}}\end{array}\right)+t^8\left(\begin{array}{c}-{\displaystyle \frac{1}{315}}\\ {\displaystyle \frac{1}{315}}\end{array}\right)+t^{10}\left(\begin{array}{c}{\displaystyle \frac{2}{14175}}\\ -{\displaystyle \frac{2}{14175}}\end{array}\right).$$

(108)

One can easily check that the above approximation forms the first terms of the Maclaurin power series of the exact solution $\mathbf{x}\left(t\right)={({\sin }^2\left(t\right),{\cos }^2\left(t\right))}^T$. It is worth noting here that we obtain the exact values of the successive DTM expansion coefficients. It is known that by increasing the number of terms in the series solution, the error decreases, and the convergence interval of the truncated series increases (Figure 7 and Figure 8).

Example 6.

In this example, we consider the following fully implicit nonlinear system of second-order ordinary differential equations:

$$\begin{array}{ccc}\hfill {\mathbf{x}}_1^{″}+3{\mathbf{x}}_2^{″}-{\displaystyle \frac{1}{1+{({\mathbf{x}}_1^{″}+{\mathbf{x}}_2^{″})}^2}}+2({\mathbf{x}}_2-1)\sin \left(t\right)-2{\mathbf{x}}_2^{\prime }\cos \left(t\right)+1& =& 0,\hfill \end{array}$$

(109)

$$\begin{array}{ccc}\hfill 3{\mathbf{x}}_1^{″}+{\mathbf{x}}_2^{″}+{\displaystyle \frac{1}{1+{({\mathbf{x}}_1^{″}+{\mathbf{x}}_2^{″})}^2}}+2({\mathbf{x}}_1-1)\sin \left(t\right)-2{\mathbf{x}}_1^{\prime }\cos \left(t\right)-1& =& 0,\hfill \end{array}$$

(110)

where $t>0$. The sought solution is $\mathbf{x}={\left({\mathbf{x}}_1,{\mathbf{x}}_2\right)}^T$. The differential system (109) and (110) has the form (1), with $n=2$, and $m=2$. This system is subject to the following initial conditions:

$$\mathbf{x}\left(0\right)=\left(\begin{array}{c}2\\ 0\end{array}\right),{\mathbf{x}}^{\prime }\left(0\right)=\left(\begin{array}{c}1\\ -1\end{array}\right).$$

(111)

We should emphasize here that Maple 2024 software cannot solve this system. However, our technique has successfully obtained the exact solution of this implicit system in power series. By applying Algorithm 2 to (109) and (110) with initial condition (111), we can recursively determine the DTM expansion coefficients $x_k$ for $k\ge 3$. According to the initial condition (111), we have the following:

$$x_0=\left(\begin{array}{c}2\\ 0\end{array}\right),x_1=\left(\begin{array}{c}1\\ -1\end{array}\right).$$

(112)

Then, we determine the differential transform $x_2={\left(x_{21},x_{22}\right)}^T$ from the following algebraic system:

$$\mathbf{F}(2x_2,x_1,x_0,0)=\mathbf{0}.$$

(113)

The above system has a unique solution $x_2$ and simplifies to the following:

$$\begin{array}{ccc}\hfill 2x_{21}+6x_{22}-{\displaystyle \frac{1}{1+{(x_{21}+x_{22})}^2}}+3& =& 0,\hfill \end{array}$$

(114)

$$\begin{array}{ccc}\hfill 6x_{21}+2x_{22}+{\displaystyle \frac{1}{1+{(x_{21}+x_{22})}^2}}-3& =& 0.\hfill \end{array}$$

(115)

Solving this algebraic system, we obtain the following:

$$x_2=\left(\begin{array}{c}{x}_{21}\\ x_{22}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{2}}\\ -{\displaystyle \frac{1}{2}}\end{array}\right).$$

(116)

Then, we compute the Jacobian $J:={\mathbf{F}}_{{\mathbf{x}}^{″}}$ of the function $\mathbf{F}$ with respect to ${\mathbf{x}}^{″}$ at the point $(x_2,x_1,x_0,0)$. We obtain the following:

$$J=\left(\begin{array}{cc}1& 3\\ 3& 1\end{array}\right).$$

(117)

Using Theorem 2, we have the following recursion for the differential transform: $x_k$

$$x_k=-{\left(k(k-1)J\right)}^{-1}{r}_{k-2},k\ge 3,$$

(118)

where

$$r_{k-2}:=F_{k-2}(2x_2,\dots ,(k-2)(k-1)x_{k-1},0;2x_2,\dots ,(k-1)x_{k-1};x_1,\dots ,\dots ,x_{k-2}).$$

(119)

Using Algorithm 2, we compute the inverse of the matrix and the right-hand sides $r_{k-2}$ for $k\ge 3$ and deduce the following:

For $k=3$, we have the following:

$$x_3=-{\left(6J\right)}^{-1}{r}_1=\left(\begin{array}{cc}{\displaystyle \frac{1}{48}}& -{\displaystyle \frac{1}{16}}\\ -{\displaystyle \frac{1}{16}}& {\displaystyle \frac{1}{48}}\end{array}\right)\left(\begin{array}{c}0\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(120)

For $k=4$, we have the following:

$$x_4=-{\left(12J\right)}^{-1}{r}_2=\left(\begin{array}{cc}{\displaystyle \frac{1}{96}}& -{\displaystyle \frac{1}{32}}\\ -{\displaystyle \frac{1}{32}}& {\displaystyle \frac{1}{96}}\end{array}\right)\left(\begin{array}{c}-3\\ 3\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{8}}\\ {\displaystyle \frac{1}{8}}\end{array}\right).$$

(121)

For $k=5$, we have the following:

$$x_5=-{\left(20J\right)}^{-1}{r}_3=\left(\begin{array}{cc}{\displaystyle \frac{1}{160}}& -{\displaystyle \frac{3}{160}}\\ -{\displaystyle \frac{3}{160}}& {\displaystyle \frac{1}{160}}\end{array}\right)\left(\begin{array}{c}-{\displaystyle \frac{8}{3}}\\ {\displaystyle \frac{8}{3}}\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{15}}\\ {\displaystyle \frac{1}{15}}\end{array}\right).$$

(122)

For $k=6$, we have the following:

$$x_6=-{\left(30J\right)}^{-1}{r}_4=\left(\begin{array}{cc}{\displaystyle \frac{1}{240}}& -{\displaystyle \frac{1}{80}}\\ -{\displaystyle \frac{1}{80}}& {\displaystyle \frac{1}{240}}\end{array}\right)\left(\begin{array}{c}-{\displaystyle \frac{1}{4}}\\ {\displaystyle \frac{1}{4}}\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{1}{240}}\\ {\displaystyle \frac{1}{240}}\end{array}\right).$$

(123)

For $k=7$, we have the following:

$$x_7=-{\left(42J\right)}^{-1}{r}_5=\left(\begin{array}{cc}{\displaystyle \frac{1}{336}}& -{\displaystyle \frac{1}{112}}\\ -{\displaystyle \frac{1}{112}}& {\displaystyle \frac{1}{336}}\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{14}{15}}\\ -{\displaystyle \frac{14}{15}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{90}}\\ -{\displaystyle \frac{1}{90}}\end{array}\right).$$

(124)

For $k=8$, we have the following:

$$x_8=-{\left(56J\right)}^{-1}{r}_6=\left(\begin{array}{cc}{\displaystyle \frac{1}{448}}& -{\displaystyle \frac{3}{448}}\\ -{\displaystyle \frac{3}{448}}& {\displaystyle \frac{1}{448}}\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{217}{360}}\\ -{\displaystyle \frac{217}{360}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{31}{5760}}\\ -{\displaystyle \frac{31}{5760}}\end{array}\right).$$

(125)

For $k=9$, we have the following:

$$x_9=-{\left(72J\right)}^{-1}{r}_7=\left(\begin{array}{cc}{\displaystyle \frac{1}{576}}& -{\displaystyle \frac{1}{192}}\\ -{\displaystyle \frac{1}{192}}& {\displaystyle \frac{1}{576}}\end{array}\right)\left(\begin{array}{c}{\displaystyle \frac{8}{315}}\\ -{\displaystyle \frac{8}{315}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{1}{5670}}\\ -{\displaystyle \frac{1}{5670}}\end{array}\right).$$

(126)

For $k=10$, we have the following:

$$x_{10}=-{\left(90J\right)}^{-1}{r}_8=\left(\begin{array}{cc}{\displaystyle \frac{1}{720}}& -{\displaystyle \frac{1}{240}}\\ -{\displaystyle \frac{1}{240}}& {\displaystyle \frac{1}{720}}\end{array}\right)\left(\begin{array}{c}-{\displaystyle \frac{2951}{20160}}\\ {\displaystyle \frac{2951}{20160}}\end{array}\right)=\left(\begin{array}{c}-{\displaystyle \frac{2951}{3628800}}\\ {\displaystyle \frac{2951}{3628800}}\end{array}\right).$$

(127)

Finally, we construct the following approximate solution with $K=10$:

$$\begin{array}{ccc}\hfill \mathbf{x}\left(t\right)& =& \left(\begin{array}{c}2\\ 0\end{array}\right)+t\left(\begin{array}{c}1\\ -1\end{array}\right)+t^2\left(\begin{array}{c}{\displaystyle \frac{1}{2}}\\ -{\displaystyle \frac{1}{2}}\end{array}\right)+t^4\left(\begin{array}{c}-{\displaystyle \frac{1}{8}}\\ {\displaystyle \frac{1}{8}}\end{array}\right)+t^5\left(\begin{array}{c}-{\displaystyle \frac{1}{5}}\\ {\displaystyle \frac{1}{5}}\end{array}\right)+t^6\left(\begin{array}{c}-{\displaystyle \frac{1}{240}}\\ {\displaystyle \frac{1}{240}}\end{array}\right)+t^7\left(\begin{array}{c}{\displaystyle \frac{1}{90}}\\ -{\displaystyle \frac{1}{90}}\end{array}\right)\hfill \\ & & +t^8\left(\begin{array}{c}{\displaystyle \frac{31}{5760}}\\ -{\displaystyle \frac{31}{5760}}\end{array}\right)+t^9\left(\begin{array}{c}{\displaystyle \frac{1}{5670}}\\ -{\displaystyle \frac{1}{5670}}\end{array}\right)+t^{10}\left(\begin{array}{c}-{\displaystyle \frac{2951}{3628800}}\\ -{\displaystyle \frac{2951}{3628800}}\end{array}\right).\hfill \end{array}$$

One can easily check that the above approximation forms the first terms of the Maclaurin power series of the exact solution $\mathbf{x}\left(t\right)={(1+e^{\sin \left(t\right)},1-e^{\sin \left(t\right)})}^T$. It is worth noting here that we obtain the exact values of the successive DTM expansion coefficients. It is known that by increasing the number of terms in the series solution, the error decreases, and the convergence interval of the truncated series increases (Figure 9 and Figure 10).

5. Discussion

In this work, we introduced a multi-stage differential transform method (MsDTM) combined with Adomian polynomials as an effective semi-analytical tool for solving fully implicit nonlinear ordinary differential systems. We demonstrated the validity of two theorems that underpin our numerical algorithms. The method was successfully applied to solve implicit ordinary differential systems of both the first- and second-order, with the potential for generalization to higher-order systems. The first unknown recursion from an algebraic system is initially solved uniquely, so the solution process starts here. This algebraic system can be nonlinear; if necessary, Newton’s method can be used to compute an accurate solution. Subsequently, the differential transform method (DTM) is used to turn the implicit ordinary differential system into a linear algebraic recursion system with a constant nonsingular coefficient matrix. This linear system can then be solved using a variety of techniques. Section 3 describes the DTM algorithm, which enables us to solve the differential system in the form of a power series. The DTM technique will produce the exact solution in power series form if the differential system has an exact solution and all computations have been performed exactly. We thereby extend the interval of convergence by using the multi-stage DTM technique described in Section 3, even though the power series solution could only converge over a short interval. One of the MsDTM method’s strengths is the extension of the convergence interval. We coded the DTM and MsDTM processes using the Maple 2024 package to apply the approach. These techniques were successfully applied to six numerical examples of fully implicit nonlinear differential systems. Maple’s built-in procedures cannot resolve these problems numerically or in power series. For each of the six cases, the DTM converted the differential system into an algebraic recursion system for the power series solution’s coefficients. We used the multi-stage MsDTM to improve the convergence of the DTM solution, hence strengthening the solution that emerged from the truncated power series. One of the main advantages of our approach is that it avoids needless computations, including first- or lower-order system reduction in the differential system. This considerably reduces the computing load, increases efficiency, and simplifies the overall technique. Furthermore, the DTM computes its coefficients by a straightforward procedure that can be easily coded in computer algebra systems like Maple, Mathematica, or Matlab. Both the DTM and the MsDTM approaches handled the numerical cases solved here with effectiveness. We used the absolute error between the exact solutions and the numerical solutions to assess the accuracy of the results. Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 show comparisons between the MsDTM approximations and the exact solutions for all numerical examples. For all approximations, the absolute error was very low. If required, the error can be further minimized by splitting the solution interval into additional intervals or by raising the order of the DTM approximations. For the numerical examples under consideration, the absolute error varied between ${10}^{-11}$ and ${10}^{-8}$, the interval of solutions was $[0,15]$, the chosen parameter was $K=7$, and the average number of interval subdivisions was $N=400$. On the other hand, techniques like Runge–Kutta (RK) require that the differential system be in an explicit form. Furthermore, using implicit numerical integration methods for ODEs typically necessitates solving nonlinear algebraic problems at each iteration step, a process that can be computationally intensive. However, by avoiding these complications, the MsDTM technique offers a more effective and useful way to solve IODEs.

A comprehensive summary of the six case studies used to assess the performance of the proposed MsDTM is given in

Table 1. Summary of main results and parameters of the six case studies. MAE stands for maximum absolute error.

Case	Order	Equations	K	N	MAE	Interval
1	1	2	7	400	<$1\times {10}^{-11}$	[0, 15]
2	1	2	7	400	<$5\times {10}^{-11}$	[0, 15]
3	1	2	7	400	<$2\times {10}^{-9}$	[0, 15]
4	2	1	5	–	Exact solution	[0, Convergent]
5	2	2	10	500	<$7.5\times {10}^{-9}$	[0, 15]
6	2	2	10	400	<$2.5\times {10}^{-8}$	[0, 15]

. Using an approximation order of $K=7$ and $N=400$ subdivisions consistently yielded maximum absolute errors (MAE) $<5\times {10}^{-11}$, indicating remarkable precision over the interval [0, 15]. The results show that this was true for first-order systems (Cases 1–3). A validation benchmark, Case 4, a second-order system with an exact solution, confirms the theoretical robustness of the DTM method. Maintaining MAE below $7.5\times {10}^{-9}$ and $2.5\times {10}^{-8}$ for more complex second-order systems (Cases 5–6) required a slightly greater number of subdivisions ($N=500$ and $N=400$ respectively) and an increased approximation order of $K=10$. These results show the method’s high accuracy for solving implicit ordinary differential equations of both the first- and second-order. Furthermore, the constant performance throughout several examples emphasizes the scalability and adaptability of MsDTM, therefore enabling MsDTM to be a powerful method for solving a broad spectrum of nonlinear IODEs.

6. Conclusions

A novel semi-analytical method, combining Adomian polynomials (AP) with the multi-stage differential transform method (MsDTM), was developed for solving implicit nonlinear ODEs. Direct application of the DTM to implicit systems simplifies the solution procedure and reduces computational costs, thus eliminating the need for explicit or quasi-linear transformations. Because AP methodically expanded nonlinear terms as power series, complex nonlinearities were efficiently handled, preserving the system’s implicit nature. The MsDTM improved upon the DTM by dividing the solution domain into smaller sub-intervals and iteratively producing solutions using truncated series as initial conditions, thus achieving high accuracy over larger domains. Numerical tests demonstrated that our method is highly accurate and reliable for solving challenging first- and second-order implicit nonlinear ODEs, even those that pose difficulties for standard methods such as Runge–Kutta, producing results very close to the exact solutions. Using AP, the MsDTM provides a robust and computationally efficient method to solve fully implicit nonlinear ODEs, improving accuracy, convergence, and applicability to mathematical problems. In research and education, the application of the MsDTM can be extended to Maple, Mathematica, and Matlab.

Future work topics will include the following:

• Solving implicit nonlinear ordinary differential systems with movable singularities.
• Comparing the accuracy and efficiency between our approach and other numerical and semi-analytical techniques to facilitate evaluation and enhancement.
• Examining convergence and error analysis to strengthen the approach’s theoretical foundation and applicability.
• Modifying the approach to increase its real-world applicability by solving problems that include stochastic elements, boundary conditions, integral equations, implicit fractional differential equations, and time-varying algebraic equations.