Complexity
Volume 2018 (2018), Article ID 8269541, 11 pages
https://doi.org/10.1155/2018/8269541

Research Article

A Numerical Algorithm for Solving Higher-Order Nonlinear BVPs with an Application on Fluid Flow over a Shrinking Permeable Infinite Long Cylinder

Laila Y. Al Sakkaf,¹ Qasem M. Al-Mdallal,² and U. Al Khawaja¹

¹Department of Physics, United Arab Emirates University, P.O. Box 15551, AlAin, UAE
²Department of Mathematical Sciences, United Arab Emirates University, P.O. Box 15551, AlAin, UAE

Correspondence should be addressed to Qasem M. Al-Mdallal; q.almdallal@uaeu.ac.ae

Received 5 November 2017; Accepted 12 February 2018; Published 12 March 2018

Academic Editor: Mohamed Belhaq

Copyright © 2018 Laila Y. Al Sakkaf et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We present an efficient iterative power series method for nonlinear boundary-value problems that treats the typical divergence problem and increases arbitrarily the radius of convergence. This method is based on expanding the solution around an iterative initial point. We employ this method to study the unsteady, viscous, and incompressible laminar flow and heat transfer over a shrinking permeable cylinder. More precisely, we solve the unsteady nonlinear Navier–Stokes and energy equations after reducing them to a system of nonlinear boundary-value problems of ordinary differential equations. The present method successfully captures dual solutions for both the flow and heat transfer fields and a unique solution at a specific critical unsteadiness parameter. Comparisons with previous numerical methods and an exact solution verify the validity, accuracy, and efficiency of the present method.

1. Introduction

Numerous phenomena in engineering and applied science fields are governed by nonlinear boundary-value problems (BVPs). Therefore, BVPs have received a huge attention from mathematicians, physicists, and engineers for the sake of finding and analyzing their solutions. Generally speaking, finding the analytical solutions for nonlinear BVPs is far from trivial and often is impossible. Therefore, many numerical techniques have been developed to solve such type of problems. These methods include Adomian’s decomposition method, homotopy perturbation method, variational iteration method, optimal homotopy asymptotic method, operational matrices techniques based on various orthogonal polynomials and wavelets, finite difference method, and spectral methods; the reader is referred to [1–6] and references therein.

The fluid dynamics and heat transfer of a viscous incompressible fluid flowing past stretching surfaces, such as a sheet or tube, have attracted considerable interests of many researchers because of their importance in many industrial applications such as the quality of certain products. One of the most interesting conditions for stretching surfaces problems is the velocity at the surface, where it mainly figures the characteristics of the fluid based upon two essential factors, fluid viscosity and suction parameter. A remarkable interest of several researchers concentrated on tracking the existence of dual solutions for the flow within a certain range of unsteadiness and suction parameters [7–17]. Although, the literature reveals numerous research papers discussing the flow over a stretching sheet and moving plate [18–26], there are only few studies focusing on the problem of flowing past a stretching cylinder or tube; see [8–10] and references therein.

It is established that the differential equations describing the fluid flow BVPs are highly nonlinear and demand extremely accurate numerical schemes. In this work, we show that an iterative procedure, based on successive power series expansions, provides one such high accuracy numerical scheme. Often, using power series solutions turns to be useless because the resulting solution diverges at a finite radius of convergence. The divergence is intrinsic to the nature of the solution in the sense that it persists to exist even with an infinite power series expansion. The method presented here solves this problem showing that the radius of convergence can be delayed arbitrarily to any large value. This value could, in principle, approach infinity achieving exact solutions.

Among the many different methods of solving nonlinear differential equations [27–33], the power series method is the most straight forward and efficient [34]. It has been used as a powerful numerical scheme for many problems [35–43] including chaotic systems [44–47]. Many numerical algorithms and codes have been developed based on this method [34–36, 44–48]. However, the abovementioned finiteness of radius of convergence is a serious problem that hinders the use of this method to wide class of differential equations, in particular the nonlinear ones. For instance, the nonlinear Schrödinger equation (NLSE) with cubic nonlinearity has as a solution. Using the power series method to solve this equation produces the power series of , which is valid only for .

Recognizing that a powerful numerical scheme based on this method is already established [34–36, 44–48], we nonetheless present a thorough investigation of the error associated with this method with the aim of showing how we can systemically reduce errors to infinitesimal values while having the Central Processing Unit (CPU) time within a reasonable range. We will show robustness and efficiency of the method via a highly demanding fluid flow boundary-value problem. Therefore, solving the problem of finite radius of convergence will open the door wide for applying the power series method to much larger class of differential equations, particularly the nonlinear ones.

Briefly, the present technique is based on iterative power series expansions of the solution. The domain of the independent variable, say , is divided into a number of segments each of width , where is smaller than the radius of convergence. A power series solution is obtained by expanding the solution around the left end of the first segment using the initial conditions given with the problem. Similarly, a power series solution is obtained by expanding around the start of the second segment but now using the first series to calculate the initial conditions. This is repeated times till a solution at is obtained. In the limit and the series solution becomes an exact solution. This scheme is effectively equivalent to an iterative procedure of repeated iterative calculation of the recursion relations of the power series in the first segment. Another aim of this paper is to apply this method to study the unsteady flow and heat transfer characteristics of fluid flow over a shrinking permeable infinite long cylinder. We will show that the iterative numerical scheme resulting from this method is exceeding the efficiency of typical numerical methods used. In addition, we managed to find an exact solution which enabled us to calculate accurately the error for a finite value of number of iterations . It should be noted that the present work is a part of the Master thesis [49].

The rest of the paper is organized as follows. In Section 2, a mathematical representation of our method is illustrated using a general form of nonlinear ordinary differential equation, while henceforth we call it iterative power series method. Sections 3 and 4 display the implementation of the iterative power series method on the heat and mass transfer model. Section 5 focuses on analyzing the validity of this present technique and demonstrating its efficiency by drawing comparisons with the achieved exact analytical solution and other numerical methods. In Section 6, we analyze and discuss the properties of the solutions obtained. We end with a summary of our main conclusions in Section 7.

2. General Scheme of the Iterative Power Series (IPS) Method

In this section, we give a brief description of the present method. Consider a general ordinary differential equation of the form

with initial conditions

where is the th derivative of , are real constants, and is a known function. The factor is introduced, without loss of generality, for the constants to correspond to the coefficients of the power series expansion below. At first, we divide the interval into a number of identical segments each of width . Then we expand in a power series around the beginning of each interval, namely,

where is the power series expansion around the start of the th segment and are the coefficients of the power series. Recursion relations between the coefficients are obtained upon substituting the power series solution, (3), into the differential equation, (1), which can be expressed in terms of the first coefficients corresponding to the initial conditions

where denotes the set of coefficients . The essential idea of the IPS method is to calculate the coefficients of the th power series from the th series according to (2),

which upon using (3) reads

and is simplified to

and then imposes the condition . Here, is the binomial function. The last equation is the basis for the IPS algorithm. Starting from the initial conditions for the power series of the zeroth interval, an iterative application of (7) leads to the coefficients of the th interval, namely, which give the solution at the desired point, ,

Both analytical and numerical schemes may be deduced from this algorithm. For the numerical scheme, the value of used is inserted as a number . On the other hand, leaving as a variable results in an analytical solution in terms of a power series in which is equivalent to a functional transformation on the zeroth order series; that is, the coefficients of the th series are functional transformation of the th series. In such a case the last power series for the th interval corresponds to such that functional transformations and all power series expansions of the zeroth up to th intervals will be included in the th expansion.

The coefficient of each th expansion represents the value of the solution at , which gives a discrete representation of . Therefore, in the limit , the discrete representation turns to a continuous one and thus we conjecture that the exact solution is obtained in the limits of and

3. Heat and Mass Transfer Model

In this section we will employ the present numerical technique on heat and mass transfer over a shrinking permeable cylinder described in Zaimi et al. [8] and Elnajjar et al. [10]. For completeness, we redescribe precisely the physical model. The flow is considered an unsteady, laminar, viscous, and incompressible fluid with uniform velocity and uniform temperature over a permeable shrinking circular cylinder. The cylinder is assumed to be infinitely long and the flow has constant properties. The diameter of the cylinder is assumed to be time dependent with the radius , where is a positive constant, is the constant of expansion/contraction strength, and is the time. Clearly, the cylinder’s radius is shrinking with time if is positive and stretching with time if is negative. Notice that since the flow is axisymmetric, the flow field should be a function of the radial coordinate, , and the longitudinal. The governing equations for the unsteady and incompressible fluid without body force are the continuity, momentum, and energy equations. These equations in cylindrical coordinate system, , are given by

where and are the polar coordinates in the radial and axial directions, respectively, and are the fluid velocity components in the radial and axial directions, respectively, and is the fluid temperature. The function represents the fluid pressure and the parameters , , and denote the fluid viscosity, the fluid density, and the fluid thermal diffusivity, respectively. Notice that we assumed that there is no azimuthal velocity component. The assumed boundary conditions associated with (10) for the velocity components and the temperature are given by

where is the constant surface temperature and is a positive constant.

The similarity transformations which convert (10) into nonlinear ordinary differential equations are given by [8, 10]

where and is the similarity variable given by

In addition, it should be noted that represents the dimensionless stream function and represents the normalized temperature. Applying the above similarity transformations, (10) and the boundary conditions (11) reduce to

subject to

where is the unsteadiness parameter representing the strength of contraction/expansion, is the suction parameter, and is the Prandtl number.

Our main target is to solve (14) and (15), subject to the boundary conditions (16), using the present technique in the ranges and at . In this study, we will analyze the normalized skin friction coefficient, , and the normalized heat transfer rate, .

4. IPS Method for Heat and Mass Transfer Model

The following is a detailed implementation of the IPS method used to solve (14) and (15) subject to the boundary conditions (16). It is worth mentioning that (14) includes only the variable , while (15) includes both and . Therefore, it is more convenient to find the variable from (14) and then solve (15). Furthermore, for the sake of simplicity, we render (14) to an initial-value problem; that is,

with

where must be chosen using the shooting method [10] so that the solution satisfies the boundary condition . Notice that by setting different initial values for in the shooting method, dual solutions are obtained.

As mentioned in Section 2, we start by expanding in power series around the initial point

and the derivatives, , , and can be calculated simply by differentiating this series. Substituted in (14), the coefficients, , for 3, can be found recursively in terms of the initial conditions , , and through the recursion relations. The first two recursion relations are given by

Recalculating , , and at gives

Now, , , and play the role of the initial conditions for the next series expansion, where we expand the solution and its derivatives in power series around

Resubstituting these power series expansions in the differential equation, we get the new recursion relations . The next iterative step is to calculate and its first two derivatives at which will give the initial conditions for the new power series. Repeating this process times, the general forms of the first two recursion relations of (14) are found to be

where

In summary, the IPS procedure can be reduced to the following algorithm:

where and are the recursion relations obtained from the differential equation. We have removed the superscripts that indicate the index of the iteration for convenience. The scheme is thus described simply as follows: one starts with (26) to calculate , followed by updating the initial conditions according to (27) and then using the updated values back in (26) and so on. The procedure has to be repeated times with .

Similarly, for the energy equation

subject to

where must be chosen using the shooting method [10] so that the solution satisfies the boundary condition . We expand in power series around the same initial point

and similarly for the derivatives, and . Substituted in (15), the coefficients, , for 2, are found recursively in terms of the initial conditions and through the recursion relations. Employing the IPS method, the first two recursion relations are found to be

where

5. Validation

In this section we aim at demonstrating the performance and efficiency of the present numerical scheme. Firstly, in order to obtain accurate numerical results, we have to pay attention to the selection of the numerical algorithm parameters, , , and . To achieve this target we focus on studying the following:

with

where is updated using the shooting method to satisfy the condition . Notice that several recent studies such as [10–16] reported the existence of a critical value of (named ) at which the problem has no solution (for ), only one solution (at ), and dual solutions (for ). It is worth mentioning that we succeeded in finding the explicit analytical form of the first solution for (33) subject to (34) under a condition ; that is

This exact solution will play a crucial role in proving the advantages of our numerical scheme.

The approximate solutions of the problems (33) and (34) at three different iterations , and , together with the exact solution obtained by (35) when and , are displayed in Figure 1. It is clearly seen that increasing the number of iterations in the IPS method delays the divergence point.

Figure 1: Velocity profiles, , for the first solution; exact (solid) and approximate (dashed) solutions at different with and .

To achieve an “optimal choice” of , we solve the problem with . Table 1 shows the values of up to 50 digits corresponding to the values of . It is clearly seen that the value of stabilizes at around ; hence we choose as the optimal value for the rest of the calculations in the entire paper. It should be noted that most of the used numerical schemes for such type of problems, [8–10, 17], had chosen to represent the infinity which, definitely, gave lower order of accuracy. This conclusion can be easily tested via the exact solution (35) which gives and . However, we will only show the interval up to for the rest of the coming figures.

Table 1: Progressing of values with at and .

7
1.06032130948772355259169445286450439527129433161564

8
1.02826614125492573507123370718399515004546832373446

9
1.01244316459533191876326159442255228857229463473780

10
1.00535010609275538041830857835707790239466505001516

11
1.00221279175507729615034859048196423158532888754168

12
1.00084874204304334646862685662321915732040012124886

13
1.00030312215822976659593816307972112761442861473174

14
1.00016671718702637162776598969384662018793573810246

15
1.00003031221582297665959381630797211276144286147317

16
1.00003031221582297665959381630797211276144286147317

17
1.00003031221582297665959381630797211276144286147317

The upper bound of the error in the IPS method can be estimated as follows. At each iterative step an error of order results from terminating the power series at . This error will be magnified times due to the iterative procedure. As a result, the upper bound of the error of the IPS method is

Table 2 presents the CPU time in seconds, which is machine-dependent, versus the upper bound of the error, , of the IPS method for the first solution at and , where varies from 3 to 8.

Table 2: The upper bound of the error for the first solution at and versus the CPU time at different values of .

CPU time (seconds)

3

0.0156001

4

0.0312002

5

0.0780005

6

0.156001

7

0.296402

8

0.546004

The exact solution, (35), provides a unique possibility of calculating the error of the IPS method and comparing it with that of other numerical methods. Figure 2 presents a comparison between the error of the IPS method and the explicit Runge–Kutta method of order four (ERK4) for the problem at and . The advantages of the present technique over the other one are notable.

(a)

(b)

Figure 2: Error of the IPS method (a) and ERK4 (b) for the case of and . The error is defined as the difference between the numerical solution and the exact solution (35).

A further comparison is done in Figure 3 on the normalized skin friction coefficient, , as a function of with Zaimi et al. [8] and Elnajjar et al. [10] for the case when . Excellent agreements are obtained. It should be mentioned herein that Elnajjar et al. [10] used a combination of the implicit Runge–Kutta method and the shooting method while Zaimi et al. [8] implemented the shooting method described in the book by Jaluria and Torrance [11].

Figure 3: Normalized skin friction coefficient, , as a function of for .

It is worth mentioning here that our technique successfully showed its definite capability to exceed the machine precision which is . A clear sign on this is the systematic reduction in the error in Table 2 when compared with the exact solution. Moreover, all computations are performed with at least 14 decimal digits of precision, knowing that all computations are operated using Mathematica software 10.4 and carried out on a Lenovo PC with the following specifications: model: Z470, processor: Core(TM) i5-2430M CPU @ 2.40 GHz, system type: 64-bit, and installed memory (RAM): 4.00 GB.

6. Results and Discussion

In this section, we discuss the effects of both suction parameter, , and unsteadiness parameter, , on the velocity profile, , the normalized skin friction coefficient, , the temperature profile, , and the heat transfer rate, . The numerical simulations are conducted at a fixed Prandtl number, , while the ranges considered for the other parameters are and .

Figure 4 shows the first and second solutions of the velocity profiles for with a fixed unsteadiness parameter, . It is clearly seen that the first solution for the fluid velocity inside the boundary layer region increases as increases, while the second solution shows an opposite trend. In addition, the two solutions of the velocity profile become steeper with higher magnitudes as increases. These observations emphasize the effect of increasing the suction parameter of the cylinder’s wall which is to decrease the boundary layer thickness. Consequently, increasing the suction parameter causes an increment in the normalized skin friction coefficient for the first solution and decrement in the normalized skin friction coefficient for the second solution, as clearly shown in Figure 5. These findings are consistent with the results reported by Elnajjar et al. [10] and Zaimi et al. [8].

(a)

(b)

Figure 4: Velocity profiles, , for different values of at : first solution (a) and second solution (b).

Figure 5: Normalized skin friction coefficient, , as a function of for different values of : the first solution (solid) and second solution (dotted).

Figure 6 presents the temperature profiles of the fluid flow, , at and . It is obviously noticeable that both solutions for temperature profiles admit similar behaviour, where they become wider and more relaxed as the suction parameter decreases. These kinds of behaviour inspire us to conclude that the developed thermal boundary layer and the corresponding rate of heat transfer are decreasing as increases. However, the second solution depicts more relaxed behaviour compared with the first solution. This slight difference between the first and second temperature profiles indicates that the second solution reflects higher thermal boundary layer than the first solution and, thus, a larger rate of the heat transfer as confirmed by Figure 7.

(a)

(b)

Figure 6: Temperature profiles, , for different values of at : first solution (a) and second solution (b).

Figure 7: Heat transfer rate, , as a function of for different values of : the first solution (solid) and second solution (dotted).

The variation of both the normalized skin friction coefficient, , and the heat transfer rate, , as functions of , are shown, respectively, in Figures 5 and 7 for . The results demonstrate the existence of a critical value in the -domain at which the problem has no solution for , only one solution at , and dual solutions for . Figure 5 shows that increases as increases which is due to the increase in the surface shear stress coefficient. Moreover, we observe that is decreasing with . However, Figure 7 clearly shows that increasing will definitely increase the heat transfer rate while increasing causes a decrease in the heat transfer rate.

Figure 8 displays the first and second solutions of the velocity profiles for with a fixed value of the suction parameter, . Generally speaking, the behaviour of is very similar to the case of the variable suction parameter; that is, increasing the unsteadiness parameter produces steeper behaviour in the velocity profiles for the first solution while the second solution shows an opposite trend. In agreement with the case of the variable suction parameter in Figure 4, increasing the unsteadiness parameter will then cause a reduction in the thickness of the boundary layer.

(a)

(b)

Figure 8: Velocity profiles, , for different values of at : first solution (a) and second solution (b).

Figure 9 presents the temperature profiles of the fluid flow, , at and . Clearly, the increase in the unsteadiness parameter or the suction parameter leads to the same trend.

(a)

(b)

Figure 9: Temperature profiles, , for different values of at : first solution (a) and second solution (b).

Finally, we end up our discussion with Figure 10 which presents an overview of the solution for problems (14) and (15) subject to the boundary conditions (16) in the domain for and . The straight line in this figure represents the occurrence of unique solution of the problem.

Figure 10: An overview of the solution for problems (14) and (15) subject to the boundary conditions (16) in the domain.

7. Conclusions

In this work, we presented a numerical technique for solving nonlinear BVPs based on iterative power series solutions. We have demonstrated its efficiency and accuracy through validation against the numerical ERK4. We have shown that our method excels over the ERK4 by orders of magnitude in accuracy. Moreover, the accuracy in our method is systematically controlled such that the error can be reduced to any arbitrary small value. We successfully studied the unsteady viscous flow over a contracting cylinder using the present technique. The velocity and temperature profiles of the ordinary version of Navier–Stokes equations for different suction and unsteadiness parameters are calculated. We have obtained the exact solution of the first solution of the fluid flow under a specific condition, , and have employed it to emphasize the efficiency of the present numerical technique.

The convergence analysis of the IPS method is considered for a future work. However, we strongly believe that the method is convergent. This is conjectured by the systematic reduction in error upon increasing the number of iterations or the number of terms in the power series. We believe that this technique will serve researchers in different fields working on nonlinear systems. In particular, the technique will be very useful for systems described by nonintegrable nonlinear differential equations.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to express their sincere appreciation to the United Arab Emirates University, Al Ain, UAE, for providing the financial support with the UPAR (7) 2015 and the UPAR (4) 2016 grants.

References

A.-M. Wazwaz, “A reliable modification of Adomian decomposition method,” Applied Mathematics and Computation, vol. 102, no. 1, pp. 77–86, 1999.

A.-M. Wazwaz, “Dual solutions for nonlinear boundary value problems by the variational iteration method,” International Journal of Numerical Methods for Heat & Fluid Flow, vol. 27, no. 1, pp. 210–220, 2017.

V. Marinca and N. Herisanu, “The optimal homotopy asymptotic method: Engineering applications,” The Optimal Homotopy Asymptotic Method: Engineering Applications, pp. 1–465, 2015.

Q. M. Al-Mdallal and F. M. Mahfouz, “Heat transfer from a heated non-rotating cylinder performing circular motion in a uniform stream,” International Journal of Heat and Mass Transfer, vol. 112, pp. 147–157, 2017.

Q. M. Al-Mdallal, “Numerical simulation of viscous flow past a circular cylinder subject to a circular motion,” European Journal of Mechanics - B/Fluids, vol. 49, no. part A, pp. 121–136, 2015.

Q. M. Al-Mdallal, “A numerical study of initial flow past a circular cylinder with combined streamwise and transverse oscillations,” Computers & Fluids, vol. 63, pp. 174–183, 2012.

P. D. Weidman, D. G. Kubitschek, and A. M. J. Davis, “The effect of transpiration on self-similar boundary layer flow over moving surfaces,” International Journal of Engineering Science, vol. 44, no. 11-12, pp. 730–737, 2006.

K. Zaimi, A. Ishak, and I. Pop, “Unsteady flow due to a contracting cylinder in a nanofluid using Buongiorno's model,” International Journal of Heat and Mass Transfer, vol. 68, pp. 509–513, 2014.

T. Fang, J. Zhang, and Y. Zhong, “Note on unsteady viscous flow on the outside of an expanding or contracting cylinder,” Communications in Nonlinear Science and Numerical Simulation, vol. 17, no. 8, pp. 3124–3128, 2012.

E. J. Elnajjar, Q. M. Al-Mdallal, and F. M. Allan, “Unsteady Flow and Heat Transfer Characteristics of Fluid Flow over a Shrinking Permeable Infinite Long Cylinder,” Journal of Heat Transfer, vol. 138, no. 9, Article ID 092003, 2016.

Y. Jaluria and K. E. Torrance, Comput Heat Transfer, 2003.

J. H. Merkin, “On dual solutions occurring in mixed convection in a porous medium,” Journal of Engineering Mathematics, vol. 20, no. 2, pp. 171–179, 1986.

J. Paullet and P. Weidman, “Analysis of stagnation point flow toward a stretching sheet,” International Journal of Non-Linear Mechanics, vol. 42, no. 9, pp. 1084–1091, 2007.

S. D. Harris, D. B. Ingham, and I. Pop, “Mixed convection boundary-layer flow near the stagnation point on a vertical surface in a porous medium: brinkman model with slip,” Transport in Porous Media, vol. 77, no. 2, pp. 267–285, 2009.

A. Postelnicu and I. Pop, “Falkner-Skan boundary layer flow of a power-law fluid past a stretching wedge,” Applied Mathematics and Computation, vol. 217, no. 9, pp. 4359–4368, 2011.

A. V. Roşca and I. Pop, “Flow and heat transfer over a vertical permeable stretching/shrinking sheet with a second order slip,” International Journal of Heat and Mass Transfer, vol. 60, no. 1, pp. 355–364, 2013.

A. Ishak, R. Nazar, and I. Pop, “Uniform suction/blowing effect on flow and heat transfer due to a stretching cylinder,” Applied Mathematical Modelling, vol. 32, no. 10, pp. 2059–2066, 2008.

N. Bachok, A. Ishak, and I. Pop, “Boundary layer stagnation-point flow and heat transfer over an exponentially stretching/shrinking sheet in a nanofluid,” International Journal of Heat and Mass Transfer, vol. 55, no. 25-26, pp. 8122–8128, 2012.

A. Ishak, Y. Y. Lok, and I. Pop, “Stagnation point flow over a shrinking sheet in micropolar fluid,” Chemical Engineering Communications, vol. 197, no. 11, pp. 1417–1427, 2010.

N. A. Yacob, A. Ishak, and I. Pop, “Melting heat transfer in boundary layer stagnation-point flow towards a stretching/shrinking sheet in a micropolar fluid,” Computers & Fluids, vol. 47, no. 1, pp. 16–21, 2011.

C. Y. Wang, “Stagnation flow towards a shrinking sheet,” International Journal of Non-Linear Mechanics, vol. 43, no. 5, pp. 377–382, 2008.

N. Bachok, A. Ishak, and I. Pop, “Unsteady boundary-layer flow and heat transfer of a nanofluid over a permeable stretching/shrinking sheet,” International Journal of Heat and Mass Transfer, vol. 55, no. 7-8, pp. 2102–2109, 2012.

Y. Zhong and T. Fang, “Unsteady stagnation-point flow over a plate moving along the direction of flow impingement,” International Journal of Heat and Mass Transfer, vol. 54, no. 15-16, pp. 3103–3108, 2011.

K. Bhattacharyya, S. Mukhopadhyay, and G. C. Layek, “Slip effects on boundary layer stagnation-point flow and heat transfer towards a shrinking sheet,” International Journal of Heat and Mass Transfer, vol. 54, no. 1–3, pp. 308–313, 2011.

F. Aman and A. Ishak, “Mixed convection boundary layer flow towards a vertical plate with a convective surface boundary condition,” Mathematical Problems in Engineering, vol. 2012, Article ID 453457, 11 pages, 2012.

N. Bachok, A. Ishak, and I. Pop, “Stagnation point flow toward a stretching/shrinking sheet with a convective surface boundary condition,” Journal of The Franklin Institute, vol. 350, no. 9, pp. 2736–2744, 2013.

A. D. Polyanin and V. F. Zaitsev, Handbook of Nonlinear Partial Differential Equations, CRC press, 2004.

C. S. Gardner, J. M. Greene, M. D. Kruskal, and R. M. Miura, “Method for solving the Korteweg-deVries equation,” Physical Review Letters, vol. 19, no. 19, pp. 1095–1097, 1967.

P. D. Lax, “Integrals of nonlinear equations of evolution and solitary waves,” Communications on Pure and Applied Mathematics, vol. 21, no. 5, pp. 467–490, 1968.

M. A. Salle and V. B. Matveev, Darboux transformations and solitons, 1991.

R. Hirota, “The Direct Method in Soliton Theory,” in Solitons, pp. 157–176, Springer, Berlin, Heidelberg, 1980.

G. Adomian, Solving Frontier Problems of Physics: The Decomposition Method, Springer, The Netherlands, 1994.

S. Liao and Y. Tan, “A general approach to obtain series solutions of nonlinear differential equations,” Studies in Applied Mathematics, vol. 119, no. 4, pp. 297–354, 2007.

R. Barrio, M. Rodríguez, A. Abad, and F. Blesa, “Breaking the limits: the Taylor series method,” Applied Mathematics and Computation, vol. 217, no. 20, pp. 7940–7954, 2011.

G. Corliss and Y. F. Chang, “Solving ordinary differential equations using Taylor series,” ACM Transactions on Mathematical Software, vol. 8, no. 2, pp. 114–144, 1982.

Y. F. Chang and G. Corliss, “ATOMFT: solving ODEs and DAEs using Taylor series,” Computers & Mathematics with Applications, vol. 28, no. 10-12, pp. 209–233, 1994.

J. D. Pryce, “Solving high-index DAEs by Taylor series,” Numerical Algorithms, vol. 19, no. 1-4, pp. 195–211, 1998.

R. Barrio, “Performance of the Taylor series method for ODEs/DAEs,” Applied Mathematics and Computation, vol. 163, no. 2, pp. 525–545, 2005.

N. S. Nedialkov and J. D. Pryce, “Solving differential-algebraic equations by TAYlor series. I. Computing TAYlor coefficients,” BIT Numerical Mathematics, vol. 45, no. 3, pp. 561–591, 2005.

N. S. Nedialkov and J. D. Pryce, “Solving differential-algebraic equations by Taylor series. (II): Computing the system Jacobian,” BIT Numerical Mathematics, vol. 47, no. 1, pp. 121–135, 2007.

N. S. Nedialkov and J. D. Pryce, “Solving differential algebraic equations by Taylor series. (III): THE DAETS code,” JNAIAM. Journal of Numerical Analysis, Industrial and Applied Mathematics, vol. 3, no. 1-2, pp. 61–80, 2008.

N. S. Nedialkov, K. R. Jackson, and G. F. Corliss, “Validated solutions of initial value problems for ordinary differential equations,” Applied Mathematics and Computation, vol. 105, no. 1, pp. 21–68, 1999.

W. Tucker, “A Rigorous ODE Solver and Smale's 14th Problem,” Foundations of Computational Mathematics, vol. 2, no. 1, pp. 53–117, 2002.

R. Barrio and F. Blesa, “Systematic search of symmetric periodic orbits in 2DOF Hamiltonian systems,” Chaos, Solitons & Fractals, vol. 41, no. 2, pp. 560–582, 2009.

R. Barrio, F. Blesa, and S. Serrano, “Periodic, escape and chaotic orbits in the Copenhagen and the (n+ 1) body ring problems,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 5, pp. 2229–2238, 2009.

R. Barrio, “Painting chaos: a gallery of sensitivity plots of classical problems,” International Journal of Bifurcation and Chaos, vol. 16, no. 10, pp. 2777–2798, 2006.

R. Barrio and S. Serrano, “A three-parametric study of the Lorenz model,” Physica D: Nonlinear Phenomena, vol. 229, no. 1, pp. 43–51, 2007.

À. Jorba and M. Zou, “A software package for the numerical integration of ODEs by means of high-order Taylor methods,” Experimental Mathematics, vol. 14, no. 1, pp. 99–117, 2005.

L. Y. Al Sakkaf, Iterative power series solutions of nonlinear partial differential equations with applications in nonlinear Science [M.S. thesis], United Arab Emirates University, Al Ain, United Arab Emirates, 2018.