Electronic Journal of Differential Equations, Vol. 2021 (2021), No. 58, pp. 1–20. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu AN ELEMENTARY METHOD FOR OBTAINING GENERAL SOLUTIONS TO SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS MARIANITO R. RODRIGO Abstract. An analytical method is proposed for finding the general solution of a system of ordinary differential equations (ODEs). The general solution is expressed as a series which in some cases can be summed to give an ex- pression in closed form. A sufficient condition for the series to converge is derived. Illustrative examples are given for scalar first-order ODEs (Riccati, Abel, homogeneous, Bernoulli, linear, separable) and for higher order ODEs (Airy, linear oscillator, Liénard, van der Pol). The method relies only on a calculus background. 1. Introduction Consider a system of ordinary differential equations (ODEs) dyj dx = fj(x, y1, y2, . . . , yd), j = 1, 2, . . . , d, where yj = yj(x) for j = 1, 2, . . . , d. In matrix form this system can be expressed as dy dx = f(x,y), (1.1) where y = (y1, y2, . . . , yd) and f = (f1, f2, . . . , fd). An important special case of (1.1) is the scalar first-order ODE (i.e. d = 1) dy dx = f(x, y). (1.2) It is well known that a scalar mth-order ODE can always be rewritten in the form (1.1). There are a number of techniques for finding particular solutions of (1.1), most notably using symmetry methods [2, 3, 4, 7]. A good compilation of techniques for obtaining exact and approximate solutions of ODEs is given in [9]. The reader is also referred to the encyclopedic treatises [5, 8] on exact solutions of ODEs. It is to be noted that in most cases the solution techniques tend to be ad hoc, i.e. specific to the type of ODE being considered. However, Lie group analysis can be quite general but is arguably non-elementary. 2010 Mathematics Subject Classification. 34A05, 34A25, 34A34. Key words and phrases. Ordinary differential equation; general solution; analytical method; exact solution. c©2021. This work is licensed under a CC BY 4.0 license. Submitted March 17, 2021. Published June 26, 2021. 1 2 M. R. RODRIGO EJDE-2021/58 The main objective of this article is to propose an elementary analytical method to find the general solution of (1.1). The term “elementary” here means a calculus background is sufficient. By “general solution” we mean a solution of (1.1) with d arbitrary constants. The general solution is expressed as a series which in some cases can be summed to obtain a closed-form expression. A convergence condition for the series is derived and yields an approximate solution when a closed form is difficult to obtain. A secondary objective of this article is pedagogical in nature. In a first course on differential equations, students are taught “recipes” for solving different types of ODEs. A solution technique may be valid for one ODE but not for another. This is of course to be expected. However, it would be useful for the students to realize that many of these techniques can be considered under a unified framework. The outline of this paper is as follows. The derivation of the proposed solution method is given in Section 2. Illustrative examples for the scalar ODE (1.2) are provided in Section 3, namely (i) separable equations, (ii) Bernoulli and first-order linear equations, (iii) homogeneous equations, (iv) Riccati equations and (v) Abel equations of the first and second kind. In Section 4 we apply the solution method to (i) second-order linear equations and (ii) Liénard equations. Brief concluding remarks are given in Section 5. 2. Derivation of the solution method Before giving the derivation of the method for the system of ODEs (1.1), it is worthwhile to look at the scalar ODE (1.2) to elucidate the underlying idea. The general solution of (1.2) can be expressed implicitly by z(x, y) = c, where c ∈ R. This can be viewed as a level curve of the function z = z(x, y). We wish to define sequences (zn)∞n=0 and (Sn)∞n=0 of functions zn = zn(x, y) and Sn = zn(x, y), respectively, such that Sn(x, y) = z0(x, y) + z1(x, y) + · · ·+ zn(x, y) and lim n→∞ Sn(x, y) = z(x, y) = c. Thus Sn(x, y) is an approximation to z(x, y). This is reminiscent of the idea in [6] but a different convergence criterion will be derived here that is amenable to a generalization to systems of ODEs. Returning now to (1.1), we can express the general solution as z(j)(x,y) = cj , j = 1, 2, . . . , d, where c1, c2, . . . , cd ∈ R. These can be thought of as level surfaces of the functions z(j) = z(j)(x,y) for j = 1, 2, . . . , d. Define the operator ∇y = ( ∂ ∂y1 , ∂ ∂y2 , . . . , ∂ ∂yd ) . Using (1.1), we see that ∂z(j) ∂x +∇yz (j) · f = 0, j = 1, 2, . . . , d. (2.1) EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 3 Define the sequences (z (j) n )∞n=0, where z (j) n = z (j) n (x,y) for each j = 1, 2, . . . , d, such that ∂z (j) n ∂x = −∇yz (j) n−1 · f , n ≥ 1. (2.2) Assume that z (j) 0 for j = 1, 2, . . . , d are given. Define also the sequences (S (j) n )∞n=0, where S (j) n = S (j) n (x,y) for each j = 1, 2, . . . , d, such that S(j) n (x,y) = z (j) 0 (x,y) + z (j) 1 (x,y) + · · ·+ z(j)n (x,y) = n∑ k=0 z (j) k (x,y), j = 1, 2, . . . , d. (2.3) Since S (j) n (x,y) is desired to be an approximation to z(j)(x,y), it must be true that lim n→∞ S(j) n (x,y) = z(j)(x,y), j = 1, 2, . . . , d. Therefore the exact general solution of (1.1) is given implicitly by ∞∑ k=0 z (j) k (x,y) = cj , j = 1, 2, . . . , d. (2.4) To evaluate how “close” S (j) n (x,y) is to z(j)(x,y), let us compare the gradient vectors associated with S (j) n (x,y) and z(j)(x,y) at each (x,y). Observing (2.1), the following conditions must hold: lim n→∞ [∂S(j) n ∂x (x,y) +∇yS (j) n (x,y) · f(x,y) ] = 0, j = 1, 2, . . . , d. (2.5) It can be seen from (2.3) and (2.2) that ∂S (j) n ∂x +∇yS (j) n · f = ∂z (j) 0 ∂x + n∑ k=1 ∂z (j) k ∂x +∇yz (j) 0 · f + n∑ k=1 (∇yz (j) k · f) = ∂z (j) 0 ∂x +∇yz (j) 0 · f + n∑ k=1 (∇yz (j) k · f −∇yz (j) k−1 · f). But the last term on the right-hand side is a telescoping sum that collapses to n∑ k=1 (∇yz (j) k · f −∇yz (j) k−1 · f) = ∇yz (j) n · f −∇yz (j) 0 · f . Thus ∂S (j) n ∂x +∇yS (j) n · f = ∂z (j) 0 ∂x +∇yz (j) n · f . Therefore the convergence condition (2.5) is equivalent to lim n→∞ [∂z(j)0 ∂x (x,y) +∇yz (j) n (x,y) · f(x,y) ] = 0, j = 1, 2, . . . , d. (2.6) In general, (2.6) may be valid only for all (x,y) belonging to some region in Rd+1. Remark 2.1. If z (j) 0 does not depend on x for j = 1, 2, . . . , d, then (2.6) becomes lim n→∞ ∇yz (j) n (x,y) · f(x,y) = 0, j = 1, 2, . . . , d. (2.7) 4 M. R. RODRIGO EJDE-2021/58 In particular, when d = 1, (2.7) simplifies further to lim n→∞ ∂zn ∂y (x, y) = 0. This is different from the convergence condition obtained in [6], where lim n→∞ (∂zn/∂y)(x, y) (∂Sn/∂y)(x, y) = 0. The reason is that the tangent vector was considered in [6], while here the gradient vector is used. The latter allows the solution method to be generalizable to systems of ODEs. 3. Examples: scalar case In this section several illustrative examples are given for some classes of scalar first-order ODEs (i.e. d = 1). 3.1. Separable equations. Consider the separable equation dy dx = f(x, y) = g(x) h(y) . (3.1) Take z0(x, y) = ∫ h(y) dy in (2.2), so that ∂z1 ∂x = −g(x) h(y) ∂z0 ∂y or z1(x, y) = − ∫ g(x) dx. Then (2.2) yields zn(x, y) = 0 for n ≥ 2 and therefore Sn(x, y) = z0(x, y) + z1(x, y), n ≥ 1. Equation (2.4) gives the general solution of the separable equation (3.1) as∫ h(y) dy − ∫ g(x) dx = c, as is well known. It is obvious that the convergence condition in (2.7) is always satisfied. 3.2. Bernoulli and first-order linear equations. Let us take a look at the ODE dy dx = f(x, y) = g(x)y + h(x)yp. (3.2) Without loss of generality, we may assume that p 6= 1, otherwise (3.2) reduces to a separable equation. Equatiion (3.2) is a first-order linear equation when p = 0 and a Bernoulli equation when p 6= 0. Choose z0(x, y) = y1−p. We claim that zn(x, y) = Fn(x)y1−p +Gn(x), n ≥ 1, F0(x) = 1, G0(x) = 0, (3.3) where Fn = Fn(x) and Gn = Gn(x) for n ≥ 1 are to be determined. Substituting (3.3) into (2.2), we have F ′n(x)y1−p +G′n(x) = (p− 1)g(x)Fn−1(x)y1−p + (p− 1)h(x)Fn−1(x), which implies that F ′n(x) = (p− 1)g(x)Fn−1(x), G′n(x) = (p− 1)h(x)Fn−1(x), n ≥ 1. EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 5 By induction on n, we can show that Fn(x) = (p− 1)n n! [ ∫ g(x) dx ]n , n ≥ 1 and therefore Gn(x) = ∫ (p− 1)n (n− 1)! [ ∫ g(x) dx ]n−1 h(x) dx. Equation (3.3) becomes zn(x, y) = (p− 1)n n! [ ∫ g(x) dx ]n y1−p + ∫ (p− 1)n (n− 1)! [ ∫ g(x) dx ]n−1 h(x) dx, n ≥ 1 and so (2.3) yields Sn(x, y) = y1−p + y1−p n∑ k=1 (p− 1)k k! [ ∫ g(x) dx ]k + n∑ k=1 ∫ (p− 1)k (k − 1)! [ ∫ g(x) dx ]k−1 h(x) dx = y1−p n∑ k=0 (p− 1)k k! [ ∫ g(x) dx ]k + n−1∑ k=0 (p− 1)k+1 k! [ ∫ g(x) dx ]k h(x) dx. As n→∞, we see that exp ( (p− 1) ∫ g(x) dx ) y1−p + (p− 1) ∫ exp ( (p− 1) ∫ g(x) dx ) h(x) dx = c. Solving for y above recovers the well-known solutions to the Bernoulli and first-order linear equations in one go. By comparison, the standard technique is to first solve a first-order linear ODE using an integrating factor, and then transform a Bernoulli equation to a first-order linear equation. Note that the convergence condition (2.7) holds since lim n→∞ ∂zn ∂y (x, y) = (1− p)y−p lim n→∞ 1 n! [ (p− 1) ∫ g(x) dx ]n = 0. 3.3. Homogeneous equations. A homogeneous equation has the form dy dx = f(x, y) = g (y x ) . (3.4) Take z0(x, y) = − log(|y|). Define Fn(u) = − ∫ u−n−1g(u)n du, Fn−1(u) = − ∫ u−ng(u)n−1 du, G(u) = ∫ 1 g(u)− u du. Then it is straightforward to show that F ′n(u) = −u−n−1g(u)n, F ′n−1(u) = −u−ng(u)n−1, F ′n(u) = g(u) u F ′n−1(u). (3.5) 6 M. R. RODRIGO EJDE-2021/58 For later use we note that − log(|u|) + ∞∑ k=1 Fk(u) = − ∞∑ k=0 ∫ 1 u [g(u) u ]k du = − ∫ 1 u 1 1− g(u)/u du = G(u), (3.6) where we used the geometric series expansion and assumed that |g(u)| < |u|. We claim that zn(x, y) = Fn (y x ) , n ≥ 1. Substituting this into (2.2) and recalling (3.5) gives ∂zn ∂x + g (y x )∂zn−1 ∂y = − y x2 F ′n (y x ) + 1 x g (y x ) F ′n−1 (y x ) = 0. This proves the claim. The left-hand side of (2.4) can be written with the aid of (3.6) as − log(|y|) + ∞∑ k=1 Fk (y x ) = − log ( |y x | ) + ∞∑ k=1 Fk (y x ) − log(|x|) = G (y x ) − log(|x|). Hence the general solution of the homogeneous equation (3.4) from (2.4) is G (y x ) − log(|x|) = c, G(u) = ∫ 1 g(u)− u du. (3.7) By comparison, (3.7) can be also be obtained by a dependent variable transforma- tion and solving the resulting separable equation. The convergence condition (2.7) is lim n→∞ ∂zn ∂y (x, y) = −1 y lim n→∞ [g(y/x) y/x ]n = 0, provided |g(y/x)| < |y/x|, which is what we assumed for g. Note, however, that (3.7) is valid even without this assumption. 3.4. Riccati equations. A general Riccati equation has the form (see [8, Sec- tion 1.2.1]) dy dx = g2(x)y2 + g1(x)y + g0(x), which can be transformed to a second-order linear ODE with variable coefficients. A general Riccati equation can also be transformed to the canonical form [8] dy dx = f(x, y) = y2 + g(x). (3.8) Unlike the previous examples, the solution to (3.8) for an arbitrary g is not known. Take z0(x, y) = y. We claim that z2n−1(x, y) = n∑ j=0 Fn,j(x)y2n−2j , z2n(x, y) = n∑ j=0 Gn,j(x)y2n−2j+1, n ≥ 1, (3.9) EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 7 where the functions Fn,j = Fn,j(x) and Gn,j = Gn,j(x) for n ≥ 1 and 0 ≤ j ≤ n are to be determined. Define G0,0(x, y) = 1. We infer from the second equation of (3.9) that z2n−2(x, y) = n−1∑ j=0 Gn−1,j(x)y2n−2j−1, n ≥ 2. (3.10) Equation (2.2) is ∂zm ∂x = −[y2 + g(x)] ∂zm−1 ∂y , m ≥ 1. (3.11) From (3.11) we see that z1(x, y) = −xy2 − ∫ g(x) dx. Therefore we deduce from (3.9) that F1,0(x) = −x, F1,1(x) = − ∫ g(x) dx. Moreover, (3.11) gives z2(x, y) = x2y3 + 2y ∫ xg(x) dx. Thus (3.9) yields G1,0(x) = x2, G1,1(x) = 2 ∫ xg(x) dx. Suppose that m = 2n in (3.11), i.e. ∂z2n ∂x = −[y2 + g(x)] ∂z2n−1 ∂y . (3.12) Substituting (3.9) into (3.12), we have n∑ j=0 G′n,j(x)y2n−2j+1 = −[y2 + g(x)] n∑ j=0 (2n− 2j)Fn,j(x)y2n−2j−1 = n−1∑ j=0 (−2n+ 2j)Fn,j(x)y2n−2j+1 + n−1∑ j=0 (−2n+ 2j)g(x)Fn,j(x)y2n−2j−1. However, n−1∑ j=0 (−2n+ 2j)g(x)Fn,j(x)y2n−2j−1 = n∑ j=1 (−2n+ 2j − 2)g(x)Fn,j−1(x)y2n−2j+1, so that n−1∑ j=1 G′n,j(x)y2n−2j+1 +G′n,0(x)y2n+1 +G′n,n(x)y = n−1∑ j=1 (−2n+ 2j)Fn,j(x)y2n−2j+1 − 2nFn,0(x)y2n+1 8 M. R. RODRIGO EJDE-2021/58 + n−1∑ j=1 (−2n+ 2j − 2)g(x)Fn,j−1(x)y2n−2j+1 − 2g(x)Fn,n−1(x)y. Equating coefficients of terms of like powers gives G′n,0(x) = −2nFn,0(x), n ≥ 2, G′n,j(x) = (−2n+ 2j)Fn,j(x) + (−2n+ 2j − 2)g(x)Fn,j−1(x), n ≥ 2, 1 ≤ j ≤ n− 1, G′n,n(x) = −2g(x)Fn,n−1(x), n ≥ 2. (3.13) Now suppose that m = 2n− 1 in (3.11), i.e. ∂z2n−1 ∂x = −[y2 + g(x)] ∂z2n−2 ∂y . (3.14) Substituting (3.9) and (3.10) into (3.14) gives n∑ j=0 F ′n,j(x)y2n−2j = −[y2 + g(x)] n−1∑ j=0 (2n− 2j − 1)Gn−1,j(x)y2n−2j−2 = n−1∑ j=0 (−2n+ 2j + 1)Gn−1,j(x)y2n−2j + n−1∑ j=0 (−2n+ 2j + 1)g(x)Gn−1,j(x)y2n−2j−2. However, n−1∑ j=0 (−2n+2j+1)g(x)Gn−1,j(x)y2n−2j−2 = n∑ j=1 (−2n+2j−1)g(x)Gn−1,j−1(x)y2n−2j , which implies that n−1∑ j=1 F ′n,j(x)y2n−2j + F ′n,0(x)y2n + F ′n,n(x) = n−1∑ j=1 (−2n+ 2j + 1)Gn−1,j(x)y2n−2j + (−2n+ 1)Gn−1,0(x)y2n + n−1∑ j=1 (−2n+ 2j − 1)g(x)Gn−1,j−1(x)y2n−2j − g(x)Gn−1,n−1(x). Equating coefficients of terms of like powers gives F ′n,0(x) = (−2n+ 1)Gn−1,0(x), n ≥ 2, F ′n,j(x) = (−2n+ 2j + 1)Gn−1,j(x) + (−2n+ 2j − 1)g(x)Gn−1,j−1(x), n ≥ 2, 1 ≤ j ≤ n− 1, F ′n,n(x) = −g(x)Gn−1,n−1(x), n ≥ 2. (3.15) The recursive equations in (3.13) and (3.15) together provide determining equa- tions for Fn,j(x) and Gn,j(x) for n ≥ 2 and 0 ≤ j ≤ n. If we set Fn,0(x) = −x2n−1, Gn,0(x) = x2n, n ≥ 2, EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 9 then G′n,0(x) + 2nFn,0 = 0, F ′n,0(x) + (2n− 1)Gn−1,0(x) = 0 and the first equations in (3.13) and (3.15) hold. The other equations in (3.13) and (3.15) can be solved recursively once g has been specified. In summary, we have G0,0(x, y) = 1, G1,0(x) = x2, G1,1(x) = 2 ∫ xg(x) dx, F1,0(x) = −x, F1,1(x) = − ∫ g(x) dx, Fn,0(x) = −x2n−1, Gn,0(x) = x2n, n ≥ 2, F ′n,j(x) = (−2n+ 2j + 1)Gn−1,j(x) + (−2n+ 2j − 1)g(x)Gn−1,j−1(x), n ≥ 2, 1 ≤ j ≤ n− 1, G′n,j(x) = (−2n+ 2j)Fn,j(x) + (−2n+ 2j − 2)g(x)Fn,j−1(x), n ≥ 2, 1 ≤ j ≤ n− 1, F ′n,n(x) = −g(x)Gn−1,n−1(x), G′n,n(x) = −2g(x)Fn,n−1(x), n ≥ 2. (3.16) The general solution of the canonical Riccati equation (3.8) from (2.4) is y + ∞∑ k=1 k∑ j=0 Fk,j(x)y2k−2j + ∞∑ k=1 k∑ j=0 Gk,j(x)y2k−2j+1 = y + ∞∑ k=1 k∑ j=0 [Fk,j(x) + yGk,j(x)]y2k−2j = c. (3.17) The convergence condition (2.7) for the canonical Riccati equation (3.8) is lim n→∞ n∑ j=0 [(2n− 2j)Fn,j(x)y2n−2j−1 + (2n− 2j + 1)Gn,j(x)y2n−2j ] = 0. (3.18) 3.4.1. Particular example: canonical Riccati equation with g(x) = 0. In this case (3.16) simplifies to G0,0(x, y) = 1, G1,0(x) = x2, G1,1(x) = 0, F1,0(x) = −x, F1,1(x) = 0, Fn,0(x) = −x2n−1, Gn,0(x) = x2n, n ≥ 2, F ′n,j(x) = (−2n+ 2j + 1)Gn−1,j(x), n ≥ 2, 1 ≤ j ≤ n− 1, G′n,j(x) = (−2n+ 2j)Fn,j(x), n ≥ 2, 1 ≤ j ≤ n− 1, F ′n,n(x) = 0, G′n,n(x) = 0, n ≥ 2. It follows that Fn,0(x) = −x2n−1, Gn,0(x) = x2n and Fn,j(x) = Gn,j(x) = 0 for all n ≥ 1 and 1 ≤ j ≤ n. The general solution (3.17) of dy/dx = y2 simplifies to y − ∞∑ k=1 x2k−1y2k + ∞∑ k=1 x2ky2k+1 = y + ( y − 1 x ) ∞∑ k=1 (x2y2)k = y 1 + xy = c provided |xy| < 1. Solving for y gives the explicit general solution y = c 1− cx . Note that this remains a solution even if the the condition |xy| < 1 is removed, as a direct substitution into the ODE shows. 10 M. R. RODRIGO EJDE-2021/58 3.4.2. Particular example: canonical Riccati equation with g(x) = x. In this case (3.16) becomes G0,0(x, y) = 1, G1,0(x) = x2, G1,1(x) = 2 3 x3 dx, F1,0(x) = −x, F1,1(x) = −1 2 x2, Fn,0(x) = −x2n−1, Gn,0(x) = x2n, n ≥ 2, F ′n,j(x) = (−2n+ 2j + 1)Gn−1,j(x) + (−2n+ 2j − 1)xGn−1,j−1(x), n ≥ 2, 1 ≤ j ≤ n− 1, G′n,j(x) = (−2n+ 2j)Fn,j(x) + (−2n+ 2j − 2)xFn,j−1(x), n ≥ 2, 1 ≤ j ≤ n− 1, F ′n,n(x) = −xGn−1,n−1(x), G′n,n(x) = −2xFn,n−1(x), n ≥ 2. Then Fn,j(x) = an,jx 2n+j−1, Gn,j(x) = bn,jx 2n+j , n ≥ 1, 0 ≤ j ≤ n, where an,j and bn,j are defined recursively by a1,0 = −1, a1,1 = −1 2 , b0,0 = 1, b1,0 = 1, b1,1 = 2 3 , an,0 = −1, bn,0 = 1, n ≥ 2, an,j = −2n+ 2j + 1 2n+ j − 1 bn−1,j + −2n+ 2j − 1 2n+ j − 1 bn−1,j−1, n ≥ 2, 1 ≤ j ≤ n− 1, bn,j = −2n+ 2j 2n+ j an,j + −2n+ 2j − 2 2n+ j an,j−1, n ≥ 2, 1 ≤ j ≤ n− 1, an,n = − 1 3n− 1 bn−1,n−1, bn,n = − 2 3n an,n−1, n ≥ 2. From (3.17) the general solution of dy/dx = y2 + x is y + ∞∑ k=1 k∑ j=0 (ak,jx 2k+j−1y2k−2j + bk,jx 2k+jy2k−2j+1) = c. The convergence condition (3.18) is lim n→∞ n∑ j=0 [(2n− 2j)an,jx 2n+j−1y2n−2j−1 + (2n− 2j + 1)bn,jx 2n+jy2n−2j ] = 0. 3.5. Abel equations. An Abel equation of the first kind has the form (see [8, Section 1.4.1]) dy dx = g3(x)y3 + g2(x)y2 + g1(x)y + g0(x), which can be transformed to the canonical form [8] dy dx = f(x, y) = y3 + g(x). (3.19) Similarly, a canonical Abel equation of the second kind has the form (see [8, Section 1.3.1]) dy dx = f(x, y) = 1 + g(x)y−1. (3.20) EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 11 As with the Riccati equation (3.8), the solution to an Abel equation for an arbi- trary g is not known. Here we will just consider (3.19) since (3.20) can be similarly handled. Take z0(x, y) = y. Equation (2.2) here is ∂zn ∂x = −[y3 + g(x)] ∂zn−1 ∂y , n ≥ 1. (3.21) It is easy to see that ∂z1 ∂x = −[y3 + g(x)] or z1(x, y) = −xy3 − ∫ g(x) dx and ∂z2 ∂x = −[y3 + g(x)](−3xy2) or z2(x, y) = 3 2 x2y5 + 3y2 ∫ xg(x) dx. We claim that zn(x, y) = n−1∑ j=0 Fn,j(x)y2n−3j+1 n ≥ 2, (3.22) where Fn,j = Fn,j(x) for n ≥ 2 and 0 ≤ j ≤ n − 1 are to be determined. We see that F2,0(x) = 3 2 x2, F2,1(x) = 3 ∫ xg(x) dx. Furthermore, (3.22) implies that zn−1(x, y) = n−2∑ j=0 Fn−1,j(x)y2n−3j−1, n ≥ 3. Thus, taking n ≥ 3 and substituting (3.22) into (3.21), we obtain n−1∑ j=0 F ′n,j(x)y2n−3j+1 = −[y3 + g(x)] n−2∑ j=0 (2n− 3j − 1)Fn−1,j(x)y2n−3j−2 = n−2∑ j=0 (−2n+ 3j + 1)Fn−1,j(x)y2n−3j+1 + n−2∑ j=0 (−2n+ 3j + 1)g(x)Fn−1,j(x)y2n−3j−2. However, n−2∑ j=0 (−2n+ 3j + 1)g(x)Fn−1,j(x)y2n−3j−2 = n−1∑ j=1 (−2n+ 3j − 2)g(x)Fn−1,j−1(x)y2n−3j+1. Hence n−2∑ j=1 F ′n,j(x)y2n−3j+1 + F ′n,0(x)y2n+1 + F ′n,n−1(x)y−n+4 12 M. R. RODRIGO EJDE-2021/58 = n−2∑ j=1 (−2n+ 3j + 1)Fn−1,j(x)y2n−3j+1 + (−2n+ 1)Fn−1,0(x)y2n+1 + n−2∑ j=1 (−2n+ 3j − 2)g(x)Fn−1,j−1(x)y2n−3j+1 + (n− 5)g(x)Fn−1,n−2(x)y−n+4. Equating coefficients of terms of like powers gives F ′n,0(x) = (−2n+ 1)Fn−1,0(x), n ≥ 3, F ′n,j(x) = (−2n+ 3j + 1)Fn−1,j(x) + (−2n+ 3j − 2)g(x)Fn−1,j−1(x), n ≥ 3, 1 ≤ j ≤ n− 2, F ′n,n−1(x) = (n− 5)g(x)Fn−1,n−2(x), n ≥ 3. (3.23) Recall that F2,0(x) = 3x2/2. We claim that Fn,0(x) = an,0x n, n ≥ 3, a2,0 = 3 2 , where (an,0)∞n=3 is to be determined. Substituting into the first equation in (3.23), we see that nan,0x n−1 = (−2n+ 1)an−1,0x n−1 or an,0 = ( − 2 + 1 n ) an−1,0, n ≥ 3. It follows that an,0 = 3 2 n∏ k=3 ( − 2 + 1 k ) , n ≥ 3. The expressions Fn,j(x) for n ≥ 3 and 1 ≤ j ≤ n− 1 can be computed recursively from the second and third equations in (3.23). From (2.4) we deduce that the general solution of the canonical Abel equation of the first kind (3.19) is y − xy3 − ∫ g(x) dx+ 3 2 x2y5 + 3y2 ∫ xg(x) dx + ∞∑ n=3 [3 2 xny2n+1 n∏ k=3 ( − 2 + 1 k ) + n−1∑ j=1 Fn,j(x)y2n−3j+1 ] = c. (3.24) The convergence condition (2.7) for the canonical Abel equation of the first kind given in (3.19) is lim n→∞ n−1∑ j=0 (2n− 3j + 1)Fn,j(x)y2n−3j = 0. (3.25) 3.5.1. Particular example: canonical Abel equation of the first kind with g(x) = x. First we claim that Fn,j(x) = an,jx n+j , n ≥ 3, 1 ≤ j ≤ n− 1. Substitution of this ansatz into the second equation of (3.23) yields (n+ j)an,jx n+j−1 = (−2n+3j+1)an−1,jx n−1+j +(−2n+3j−2)xan−1,j−1x n+j−2, EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 13 so that an,j = −2n+ 3j + 1 n+ j an−1,j + −2n+ 3j − 2 n+ j an−1,j−1 , n ≥ 3, 1 ≤ j ≤ n− 2. (3.26) Moreover, substituting the ansatz into the third equation of (3.23) gives (2n− 1)an,n−1x 2n−2 = (n− 5)xan−1,n−2x 2n−3 or an,n−1 = n− 5 2n− 1 an−1,n−2. Since F2,1(x) = x3, we see that a2,1 = 1. Define bn = an−1,n−2 for n ≥ 3. Hence bn+1 = n− 5 2n− 1 bn, n ≥ 3, b3 = 1 and we deduce that bn = n−1∏ k=3 k − 5 2k − 1 , n ≥ 4. Thus an,n−1 = bn+1 = n∏ k=3 k − 5 2k − 1 , n ≥ 3. We see from (3.24) that the general solution of dy/dx = y3 + x is y − xy3 − 1 2 x2 + 3 2 x2y5 + x3y2 + ∞∑ n=3 n−1∑ j=0 an,jx n+jy2n−3j+1 = c, where an,0 = 3 2 n∏ k=3 ( − 2 + 1 k ) , an,n−1 = n∏ k=3 k − 5 2k − 1 , n ≥ 3 and an,j for n ≥ 3 and 1 ≤ j ≤ n− 2 are calculated recursively from (3.26), while the convergence condition (3.25) is lim n→∞ n−1∑ j=0 (2n− 3j + 1)an,jx n+jy2n−3j = 0. 4. Examples: system case Let us now turn to some examples of systems of ODEs with d = 2. 4.1. Second-order linear ODEs. Consider the second-order linear ODE d2y dx2 + g(x)y = 0. (4.1) The general solution of (4.1) is not known for an arbitrary g. If one nontrivial solution can be found, then another linearly independent solution can be obtained through the method of variation of parameters and hence the general solution can be determined. However, there is no general procedure for finding a particular solution when g is arbitrary. To express (4.1) as a system of the form (1.1), define y1 = y and y2 = dy/dx, yielding dy1 dx = f1(x, y1, y2) = y2, dy2 dx = f2(x, y1, y2) = −g(x)y1. (4.2) 14 M. R. RODRIGO EJDE-2021/58 Take z (1) 0 (x, y1, y2) = y1, z (2) 0 (x, y1, y2) = y2. We claim that z (1) 2n−1(x, y1, y2) = F (1) n (x)y2, z (1) 2n (x, y1, y2) = G(1) n (x)y1, n ≥ 1 and z (2) 2n−1(x, y1, y2) = F (2) n (x)y1, z (2) 2n (x, y1, y2) = G(2) n (x)y2, n ≥ 1, where F (j) n = F (j) n (x) andG (j) n = G (j) n (x) for j = 1, 2 and n ≥ 1 are to be determined such that G (1) 0 (x) = 1 and G (2) 0 (x) = 1. Equation (2.2) becomes ∂z (1) 2n−1 ∂x = −y2 ∂z (1) 2n−2 ∂y1 + g(x)y1 ∂z (1) 2n−2 ∂y2 , ∂z (1) 2n ∂x = −y2 ∂z (1) 2n−1 ∂y1 + g(x)y1 ∂z (1) 2n−1 ∂y2 , ∂z (2) 2n−1 ∂x = −y2 ∂z (2) 2n−2 ∂y1 + g(x)y1 ∂z (2) 2n−2 ∂y2 , ∂z (2) 2n ∂x = −y2 ∂z (2) 2n−1 ∂y1 + g(x)y1 ∂z (2) 2n−1 ∂y2 . It follows that ∂z (1) 1 ∂x = −y2 or z (1) 1 (x, y1, y2) = −xy2 and ∂z (1) 2 ∂x = −xg(x)y1 or z (1) 2 (x, y1, y2) = −y1 ∫ xg(x) dx. This allows us to identify F (1) 1 (x) = −x, G (1) 1 (x) = − ∫ xg(x) dx. Similarly, we see that ∂z (2) 1 ∂x = g(x)y1 or z (2) 1 (x, y1, y2) = y1 ∫ g(x) dx and ∂z (2) 2 ∂x = −y2 ∫ g(x) dx or z (2) 2 (x, y1, y2) = −y2 ∫ [ ∫ g(x) dx ] dx, from which we identify F (2) 1 (x) = ∫ g(x) dx, G (2) 1 (x) = − ∫ [ ∫ g(x) dx ] dx. More generally, for n ≥ 2, we obtain d dx (F (1) n )y2 = −y2G(1) n−1, d dx (G(1) n )y1 = g(x)y1F (1) n , d dx (F (2) n )y1 = g(x)y1G (2) n−1, d dx (G(2) n )y2 = −y2F (2) n . EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 15 These are equivalent to d2 dx2 F (1) n + g(x)F (1) n−1 = 0, n ≥ 2, F (1) 1 (x) = −x, d2 dx2 G(2) n + g(x)G (2) n−1 = 0, n ≥ 2, G (2) 1 (x) = − ∫ [ ∫ g(x) dx ] dx, G(1) n (x) = ∫ g(x)F (1) n (x) dx, F (2) n (x) = ∫ g(x)G (2) n−1(x) dx, n ≥ 1. (4.3) Therefore (2.4) in this case is y1 + y2 ∞∑ n=1 F (1) n (x) + y1 ∞∑ n=1 G(1) n (x) = y1 ∞∑ n=0 G(1) n (x) + y2 ∞∑ n=1 F (1) n (x) = c1 , y2 + y1 ∞∑ n=1 F (2) n (x) + y2 ∞∑ n=1 G(2) n (x) = y1 ∞∑ n=1 F (2) n (x) + y2 ∞∑ n=0 G(2) n (x) = c2. (4.4) But (4.4) is a linear system in y1 and y2, and whose solution is y1 = c1 ∑∞ n=0G (2) n (x)− c2 ∑∞ n=1 F (1) n (x)∑∞ n=0G (1) n (x) ∑∞ n=0G (2) n (x)− ∑∞ n=1 F (1) n (x) ∑∞ n=1 F (2) n (x) , y2 = −c1 ∑∞ n=1 F (2) n (x) + c2 ∑∞ n=0G (1) n (x)∑∞ n=0G (1) n (x) ∑∞ n=0G (2) n (x)− ∑∞ n=1 F (1) n (x) ∑∞ n=1 F (2) n (x) . Recalling that y1 = y, the general solution of the second-order linear equation (4.1) is therefore y = c1 ∑∞ n=0G (2) n (x)− c2 ∑∞ n=1 F (1) n (x)∑∞ n=0G (1) n (x) ∑∞ n=0G (2) n (x)− ∑∞ n=1 F (1) n (x) ∑∞ n=1 F (2) n (x) , (4.5) where F (j) n (x) and G (j) n (x) for j = 1, 2 and n ≥ 1 are as given recursively in (4.3). Note that G (1) 0 (x) = G (2) 0 (x) = 1 by assumption. It is easy to see that ∂z (1) 2n ∂y1 (x, y1, y2) = G(1) n (x), ∂z (1) 2n−1 ∂y1 (x, y1, y2) = 0, ∂z (2) 2n ∂y2 (x, y1, y2) = G(2) n (x), ∂z (2) 2n−1 ∂y2 (x, y1, y2) = 0. Therefore sufficient conditions for the convergence condition (2.7) to be valid are lim n→∞ G(1) n (x) = 0, lim n→∞ G(2) n (x) = 0. (4.6) 4.1.1. Particular example: second-order linear equation with g(x) = 1. Suppose that we take g(x) = 1 in (4.1), giving the classical oscillator equation [9]. It is straightforward to show that the solution of (4.3) is F (1) n (x) = (−1)n x2n−1 (2n− 1)! , G(2) n (x) = (−1)n x2n (2n)! , n ≥ 1 and G(1) n (x) = (−1)n x2n (2n)! = G(2) n (x), 16 M. R. RODRIGO EJDE-2021/58 F (2) n (x) = (−1)n−1 x2n−1 (2n− 1)! = −F (1) n (x), n ≥ 1. Moreover, we obtain ∞∑ n=0 G(1) n (x) = ∞∑ n=0 G(2) n (x) = ∞∑ n=0 (−1)n x2n (2n)! = cos(x) and ∞∑ n=1 F (1) n (x) = − ∞∑ n=1 F (2) n (x) = ∞∑ n=1 (−1)n x2n−1 (2n− 1)! = − ∞∑ n=0 (−1)n x2n+1 (2n+ 1)! = − sin(x). Using (4.5), we deduce that the general solution of d2y/dx2 + y = 0 is y = c1 cos(x) + c2 sin(x) as to be expected. Note that (4.6) is satisfied since lim n→∞ (−1)n x2n (2n)! = 0. 4.1.2. Particular example: second-order linear equation with g(x) = −x. Assume that g(x) = −x in (4.1), which yields the Airy ODE [1, 9] whose general solution is y = c1 Ai(x) + c2 Bi(x), where Ai and Bi are the Airy functions. Define the sequences (an)∞n=1 and (bn)∞n=1 by an = 1 (3n− 2)(3n− 3) an−1, n ≥ 2, a1 = −1, bn = 1 3n(3n− 1) bn−1, n ≥ 2, b1 = 1 6 . We can verify by direct substitution that the respective solutions of these difference equations are an = − n−2∏ j=0 1 (3j + 3)(3j + 4) , bn = 1 6 n−2∏ j=0 1 (3j + 5)(3j + 6) , n ≥ 2. Again by direct substitution we can show that the solution of (4.3) is F (1) n (x) = anx 3n−2, G(2) n (x) = bnx 3n, n ≥ 1 and G(1) n (x) = − 1 3n anx 3n, F (2) n (x) = − 1 3n− 1 bn−1x 3n−1, n ≥ 1. Thus from (4.5) the general solution of the Airy ODE d2y/dx2 − xy = 0 can also be expressed as y = c1R1(x) + c2R2(x), where R1(x) = ∑∞ n=0G (2) n (x)∑∞ n=0G (1) n (x) ∑∞ n=0G (2) n (x)− ∑∞ n=1 F (1) n (x) ∑∞ n=1 F (2) n (x) , EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 17 R2(x) = − ∑∞ n=1 F (1) n (x)∑∞ n=0G (1) n (x) ∑∞ n=0G (2) n (x)− ∑∞ n=1 F (1) n (x) ∑∞ n=1 F (2) n (x) . In principle, it should be possible to express Ai and Bi as linear combinations of R1 and R2, and vice versa. Equation (4.6) is satisfied if lim n→∞ [ 1 3n n−2∏ j=0 1 (3j + 3)(3j + 4) ] x3n = 0, lim n→∞ [1 6 n−2∏ j=0 1 (3j + 5)(3j + 6) ] x3n = 0. These are achievable, for instance, if −1 < x < 1. 4.2. Liénard equations. A Liénard equation [8, 9] has the form d2y dx2 + g(y) dy dx + h(y) = 0. (4.7) A particular case is the van der Pol oscillator [9], where g(y) = −α(1 − y2) for α ∈ R and h(y) = y. Exact analytical solutions to such nonlinear ODEs are not known in general. Let y1 = y and y2 = dy/dx, so that a Liénard equation has the equivalent form dy1 dx = f1(x, y1, y2) = y2, dy2 dx = f2(x, y1, y2) = −g(y1)y2 − h(y1). (4.8) Take z (1) 0 (x, y1, y2) = y1, z (2) 0 (x, y1, y2) = y2. For n ≥ 1, (2.2) becomes ∂z (1) n ∂x = −y2 ∂z (1) n−1 ∂y1 + [g(y1)y2 + h(y1)] ∂z (1) n−1 ∂y2 , ∂z (2) n ∂x = −y2 ∂z (2) n−1 ∂y1 + [g(y1)y2 + h(y1)] ∂z (2) n−1 ∂y2 . (4.9) In particular, ∂z (1) 1 ∂x = −y2, ∂z (2) 1 ∂x = g(y1)y2 + h(y1) imply that z (1) 1 (x, y1, y2) = −xy2, z (2) 1 (x, y) = x[g(y1)y2 + h(y1)], while ∂z (1) 2 ∂x = −x[g(y1)y2 + h(y1)], ∂z (2) 2 ∂x = −xy2[g′(y1)y2 + h′(y1)] + xg(y1)[g(y1)y2 + h(y1)] imply that z (1) 2 (x, y1, y2) = −1 2 x2[g(y1)y2 + h(y1)], z (2) 2 (x, y1, y2) = 1 2 x2[−g′(y1)y22 − h′(y1)y2 + g(y1)2y2 + g(y1)h(y1)]. 18 M. R. RODRIGO EJDE-2021/58 Define F0(y1, y2) = y1 and G0(y1, y2) = y2. We claim that z(1)n (x, y1, y2) = 1 n! Fn(y1, y2)xn, z(2)n (x, y1, y2) = 1 n! Gn(y1, y2)xn, n ≥ 1, (4.10) where Fn = Fn(y1, y2) and Gn = Gn(y1, y2) for n ≥ 1 are to be determined. It is not difficult to see that F1(y1, y2) = −y2, G1(y1, y2) = g(y1)y2 + h(y1), F2(y1, y2) = −g(y1)y2 − h(y1), G2(y1, y2) = −g′(y1)y22 − h′(y1)y2 + g(y1)2y2 + g(y1)h(y1). Substituting (4.10) into (4.9) and equating coefficients of like powers, we obtain Fn(y1, y2) = −y2 ∂Fn−1 ∂y1 (y1, y2) + [g(y1)y2 + h(y1)] ∂Fn−1 ∂y2 (y1, y2), Gn(y1, y2) = −y2 ∂Gn−1 ∂y1 (y1, y2) + [g(y1)y2 + h(y1)] ∂Gn−1 ∂y2 (y1, y2). (4.11) Suppose that Fn(y1, y2) = n∑ j=0 fn,j(y1)yj2, Gn(y1, y2) = n∑ j=0 gn,j(y1)yj2, n ≥ 1, (4.12) where fn,j = fn,j(y1) and gn,j = gn,j(y1) for n ≥ 1 and 0 ≤ j ≤ n are to be determined. Then f1,0(y1) = 0, f1,1(y1) = −1, g1,0(y1) = h(y1), g1,1(y1) = g(y1), f2,0(y1) = −h(y1), f2,1(y1) = −g(y1), f2,2(y1) = 0, g2,0(y1) = g(y1)h(y1), g2,1(y1) = −h′(y1) + g(y1)2, g2,2(y1) = −g′(y1). Substitution of the first ansatz in (4.12) into (4.11) gives n∑ j=0 fn,j(y1)yj2 = −y2 n−1∑ j=0 f ′n−1,j(y1)yj2 + [g(y1)y2 + h(y1)] n−1∑ j=0 jfn−1,j(y1)yj−12 = − n−1∑ j=0 f ′n−1,j(y1)yj+1 2 + n−1∑ j=1 jg(y1)fn−1,j(y1)yj2 + n−1∑ j=1 jh(y1)fn−1,j(y1)yj−12 . However, n−1∑ j=0 f ′n−1,j(y1)yj+1 2 = n∑ j=1 f ′n−1,j−1(y1)yj2 and n−1∑ j=1 jh(y1)fn−1,j(y1)yj−12 = n−2∑ j=0 (j + 1)h(y1)fn−1,j+1(y1)yj2, EJDE-2021/58 GENERAL SOLUTION TO SYSTEMS OF ODES 19 so that n−2∑ j=1 fn,j(y1)yj2 + fn,0(y1) + fn,n−1(y1)yn−12 + fn,n(y1)yn2 = − n−2∑ j=1 f ′n−1,j−1(y1)yj2 − f ′n−1,n−2(y1)yn−12 − f ′n−1,n−1(y1)yn2 + n−2∑ j=1 jg(y1)fn−1,j(y1)yj2 + (n− 1)g(y1)fn−1,n−1(y1)yn−12 + n−2∑ j=1 (j + 1)h(y1)fn−1,j+1(y1)yj2 + h(y1)fn−1,1(y1). Equating coefficients of like powers, we obtain fn,0(y1) = h(y1)fn−1,1(y1), n ≥ 3, fn,j(y1) = −f ′n−1,j−1(y1) + jg(y1)fn−1,j(y1) + (j + 1)h(y1)fn−1,j+1(y1), n ≥ 3, 1 ≤ j ≤ n− 2, fn,n−1(y1) = −f ′n−1,n−2(y1) + (n− 1)g(y1)fn−1,n−1(y1), n ≥ 3, fn,n(y1) = −f ′n−1,n−1(y1), n ≥ 3. In a similar manner, substitution of the second ansatz in (4.12) into (4.11) and replacing all occurrences of fn,j by gn,j gives gn,0(y1) = h(y1)gn−1,1(y1), gn,j(y1) = −g′n−1,j−1(y1) + jg(y1)gn−1,j(y1) + (j + 1)h(y1)gn−1,j+1(y1), n ≥ 3, 1 ≤ j ≤ n− 2, gn,n−1(y1) = −g′n−1,n−2(y1) + (n− 1)g(y1)gn−1,n−1(y1), n ≥ 3, gn,n(y1) = −g′n−1,n−1(y1), n ≥ 3. Hence from (2.4) we get the general solution of (4.8), namely y1 − xy2 + ∞∑ k=2 xk k! k∑ j=0 yj2fk,j(y1) = c1, y2 + x[g(y1)y2 + h(y1)] + ∞∑ k=2 xk k! k∑ j=0 yj2gk,j(y1) = c2. Identifying y1 = y therefore gives the general solution of the Liénard equation (4.7) from the above pair of equations. Since ∂z (1) n ∂y1 (x, y1, y2) = xn n! n∑ j=0 yj2f ′ n,j(y1), ∂z (2) n ∂y2 (x, y1, y2) = xn n! n∑ j=0 jyj−12 gn,j(y1), the convergence condition (2.7) is satisfied if lim n→∞ xn n! n∑ j=0 yj2f ′ n,j(y1) = 0, lim n→∞ xn n! n∑ j=1 jyj−12 gn,j(y1) = 0. 20 M. R. RODRIGO EJDE-2021/58 5. Concluding remarks In this article we proposed an elementary analytical method for finding the general solution of a system of ODEs given by (1.1). The exact analytical solu- tion (2.4) is expressed as a series and a convergence criterion (2.6) was derived for the series solution. The criterion determines a region in (x,y)-space for which the series will converge. By choosing different expressions for the “initial seeds” z (j) 0 for j = 1, 2, . . . , d, other regions of convergence can be obtained. In some cases, the series can be summed and a closed-form expression for the general solution can be deduced. When this is not possible, approximate analytical formulas for the general solutions can nevertheless be written down. A careful error analysis of such analytical approximations is outside the scope of this article but will be investigated in future work. We gave several illustrative examples in one and two dimensions. In particular, for Riccati, Abel and Liénard equations, analytical solutions are not known in the general case (cf. [8] for a compendium of solvable special cases). The proposed method here is an alternative to Lie symmetry analysis and is a useful addition to the applied mathematician’s toolbox. It was also shown that many known solvable ODEs (e.g. separable, first-order linear, homogeneous, Bernoulli, Airy and other second-order linear equations) can be handled with the same methodology and can thus be introduced in a differential equations class at a more elementary level. Another future research direction will be to extend the proposed method here to fractional-order ODEs and certain classes of partial differential equations. References [1] Abramowitz, M.; Stegun, I.A.; Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables, 9th edn. New York: Dover, 1972 [2] Bluman, G. W.; Anco, S.C.; Symmetry and Integration Methods for Differential Equations. Applied Mathematical Sciences vol. 154. New York: Springer, 2002 [3] Cantwell, B. J.; Introduction to Symmetry Analysis. Cambridge: Cambridge University Press, 2002 [4] Hill, J .M.; Differential Equations and Group Methods for Scientists and Engineers. Boca Raton: CRC Press, 1992. [5] Kamke, E.; Differentialgleichungen. Lösungsmethoden und Lösungen I. Gewöhnliche Differ- entialgleichungen. Leipzig: B.G. Teubner, 1971. [6] McGarvey, J. F.; Approximating the general solution of a differential equation. SIAM Review, 24 (1982) (3), pp. 333–337. [7] Olver, P. J.; Applications of Lie Groups to Differential Equations. Graduate Texts in Math- ematics vol. 107. New York: Springer, 1998. [8] Polyanin, A .D.; Zaitsev, V. F.; Handbook of Exact Solutions for Ordinary Differential Equa- tions, 2nd ed. Boca Raton: CRC Press, 2003. [9] Zwillinger, D.; Handbook of Differential Equations. New York: Academic Press, 1998. Marianito R. Rodrigo School of Mathematics and Applied Statistics, University of Wollongong, Wollon- gong, New South Wales, 2422, Australia Email address: marianito rodrigo@uow.edu.au 1. Introduction 2. Derivation of the solution method 3. Examples: scalar case 3.1. Separable equations 3.2. Bernoulli and first-order linear equations 3.3. Homogeneous equations 3.4. Riccati equations 3.5. Abel equations 4. Examples: system case 4.1. Second-order linear ODEs 4.2. Liénard equations 5. Concluding remarks References