# Preface

This section is devoted to a fascinating method for solving linear and nonlinear ordinary and partial differential equations invented by George Adomian (1922--1996). This method will be used in some other sections, so here we just give it a friendly introduction.

The Adomian decomposition method (ADM) is a systematic approximation method for solving nonlinear functional equations including ordinary and partial differential equations. It was named by Richard Bellman in honor of Adomian because it was developed from the 1970s --- 1990s by an American mathematician and aerospace engineer of Armenian descent George Adomian (1922-1996), chair of the Center for Applied Mathematics at the University of Georgia. The method is based on the assumption that the solution to a nonlinear functional equation can be represented by the infinite series $$y(x) = \sum_{n \ge 0} u_n (x) .$$ The crucial aspect of the method is the decomposition of the problem under consideration into a sequence of recursively specified auxiliary problems that can be solved explicitly. This stage utilizes the employment of so called "Adomian polynomials" to represent the nonlinear portion of the equation as a convergent series with respect to these polynomials, without actual linearization of the system. These polynomials were named after G. Adomian by Randolph Rach in 1984 (A convenient computational form for the Adomian polynomials can be found in the Journal of Mathematical Analysis and Applications, Volume 102, Issue 2, September 1984, pages 415--419). There are several known algorithms for calculating Adomian polynomials; the following publications reflect efforts made in this area:
1. Babolian, E. and Javadi, Sh., “New method for calculating Adomian polynomials”, Applied Mathematics and Computation, 2004, Volume 153, Issue 1, 25 May 2004, pages 253--259. https://doi.org/10.1016/S0096-3003(03)00629-5
2. Biazar, J., Babolian, E., Kember, G., Nouri, A., Islam, R., An alternate algorithm for computing Adomian polynomials in special cases, Applied Mathematics and Computation, 2003, Volume 138, Issues 2–3, 20 June 2003, Pages 523--529; https://doi.org/10.1016/S0096-3003(02)00174-1
3. Biazar, J. and Shafiof, S.M. A Simple Algorithm for Calculating Adomian Polynomials, International Journal of Contemporary Mathematical Sciences, 2007, Vol. 2, no. 20, pp. 975--982.
4. Chang, S.H. and Chang, I.L., A new algorithm for calculating one-dimensional differential transform of nonlinear functions, Applied Mathematics and Computation, 2008, Vol. 195, Issue 2, pp. 799--808. https://doi.org/10.1016/j.amc.2007.05.026
5. Chen, W., and Lu, Z., An algorithm for Adomian decomposition method, Applied Mathematics and Computation, 2004, Vol. 159, No. 1, (2004) pp. 221--235.
6. Choi, H.-W., and Shin, J.G., Symbolic implementation of the algorithm for calculating Adomian polynomials, Applied Mathematics and Computation, 2003, Volume 146, Issue 1, 30 December 2003, Pages 257--271; https://doi.org/10.1016/S0096-3003(02)00541-6
7. Daftardar-Gejji, V. and Jafari, H., An iterative method for solving non linear functional equations, Journal of Mathematical Analysis and Applications, 2006, vol. 316, Issue 2, pp. 753--763; https://doi.org/10.1016/j.jmaa.2005.05.009
8. Jun-Sheng Duan, Recurrence triangle for Adomian polynomials, Applied Mathematics and Computation, 2010, Vol. 216, Issue 4, pp. 1235--1241. Elsevier. 15 April 2010; https://doi.org/10.1016/j.amc.2010.02.015
9. Duan, J.-S., An efficient algorithm for the multivariable Adomian polynomials, Applied Mathematics and Computation, 2010, Vol. 217, Issue 6, pp. 2456--2467. https://doi.org/10.1016/j.amc.2010.07.046
10. Jun-Sheng Duan, Convenient analytic recurrence algorithms for the Adomian polynomials, Applied Mathematics and Computation, 2011, Vol. 217, Issue 13, pp. 6337--6348. Elsevier. 1 March 2011. https://doi.org/10.1016/j.amc.2011.01.007
11. Duan, J.-S., New recurrence algorithms for the nonclassic Adomian polynomials, Computers & Mathematics with Applications, 2011, Volume 62, Issue 8, pages 2961--2977. https://doi.org/10.1016/j.camwa.2011.07.074
12. Elsaid, A., Adomian polynomials: a powerful tool for iterative methods of series solution of nonlinear equations, Journal of Applied Analysis and Computation, Volume 2, Number 4, November 2012, pp. 381--394.
13. Fatoorehchi, H. and Abolghasemi, H. (2011b), “On calculation of Adomian polynomials by MATLAB”, Journal of Applied Computer Science and Mathematics (“Stefan cel Mare” University of Suceava), 2011, Vol. 11, No. 5, pp. 85--88.
14. Jafari, H., Ghasempoor, S., and Khaliqu, C.M., A comparison between Adomian’s polynomials and He’s polynomials for nonlinear functional equations, Mathematical Problems in Engineering, Volume 2013, Article ID 943232, 4 pages.
15. Kaliyappan, M. and Hariharan, S., Symbolic computation of Adomian polynomials based on Rach’s Rule, British Journal of Mathematics & Computer Science, Vol. 5, Issue 5, pp. 562--570, 2015, Article no.BJMCS.2015.
16. Kataria, K.K. and Vellaisamy, Simple parametrization methods for generating Adomian polynomials, Applicable Analysis and Discrete mathematics, 2016, Vol. 10, N. 01, pp. 168--185; doi: 10.2298/AADM160123001K available online at http://pefmath.etf.rs
17. Liao, S.J., On the homotopy analysis method for nonlinear problems, Applied Mathematics and Computation, 2004, Vol. 147, Issue 2, pp. 499--513; https://doi.org/10.1016/S0096-3003(02)00790-7
18. Nigam, R. A New Formulation of Adomian Polynomials, International Journal of Mathematics and Scientific Computing, Vol. 5, No 2, 2015, pp. 92--97.
19. Rach, R.C., “A new definition of the Adomian polynomials”, Kybernetes, Vol. 37 Issue: 7, 2008, pp.910--955, https://doi.org/10.1108/03684920810884342
20. Wazwaz, A.M., A new algorithm for calculating Adomian polynomials for non-linear operators, Applied Mathematics and Computation, 111 (2000), pp. 53--69.
21. Zaouagui I.N. and Badredine T., New Adomian’s Polynomials Formulas for the Non-linear and Nonautonomous Ordinary Differential Equations, Journal of Applied & Computational Mathematics, 2017, Vol 6(4), pp. doi: 10.4172/2168-9679.1000373

# Complete List of Publications using the Adomian Decomposition Method

Let me introduce some scientists (in alphabetic order) who contributed to and greatly improved the Adomian Decomposition Method to make it available to the mathematical and engineering community.

Dr. George Adomian (1922--1996) was an American mathematician, theoretical physicist, and electrical engineer of Armenian descent. He received his Ph.D. degree from UCLA. He first proposed and considerably developed the Adomian Decomposition Method (ADM) for solving nonlinear differential equations, both ordinary, and partial, deterministic and stochastic, also integral equations, algebraic and transcendental (functional), and matrix equations. He was a Distinguished Professor (academic rank), the David C. Barrow Professor of Mathematics (Chair), and the Director of the Center for Applied Mathematics at the University of Georgia, the founder and Chief Scientist of General Analytics Corporation, a winner of the 1989 Richard Bellman Prize for outstanding contributions to nonlinear stochastic analysis, and a 1988 National Academy of Sciences Scholar. He is the author of eight books and over three hundred journal papers. Dr. Adomian was the Radar Officer (naval rank of Lieutenant, the equivalent of a Captain in the Army by direct commission) aboard the USS Antietam (CV-36), and served in the Pacific Theater of Operations during the Second World War. He also attended the secret radar school at MIT that trained the first radar officers for the U.S. Navy. His obituary was prepared by Randolph Rach and appears in
R. Rach, Dr. George Adomian – distinguished scientist and mathematician, Kybernetes, Vol. 25, No. 9, 1996, 45-50.

Dr. Yves Cherruault (1937--2010) was a French mathematician. He received his Ph.D. degree from the University of Paris. He was a professor of mathematics at the University Pierre et Marie Curie, Paris, and the Director of MÉDIMAT (Laboratory of Mathematics Applied to Biomedicine). Professor Cherruault is one of the founders of the field of Biomathematics and an author of seven books and over two hundred journal papers. Dr. Cherruault developed some of the convergence theorems of the ADM.

Dr. Jun-Sheng Duan is a Chinese applied mathematician and computer scientist born in Inner Mongolia, China, in 1965. He received his Ph.D. degree from Shandong University, China. He has made extensive contributions to solutions of differential equations in mathematical physics and engineering using ADM and the Modified Decomposition Method (MDM) in collaboration with Drs. Rach and Wazwaz. He has been a professor of mathematics at the Inner Mongolia Polytechnic University, the Tianjin University of Commerce, and currently at the School of Sciences, Shanghai Institute of Technology, China. He is the author of more than eighty journal papers. Dr. Duan is an outstanding player of wei qi (the Chinese “game of go”), and ping-pong (table tennis). He enjoys photography and traveling.

Dr. Randolph Rach is a retired Army veteran and former Senior Engineer at Microwave Laboratories Inc. His experience includes research and development in microwave electronics and traveling-wave tube technology, and his research interests span nonlinear system analysis, nonlinear ordinary and partial differential equations, nonlinear integral equations and nonlinear boundary value problems. He published more than one hundred and thirty papers in applied mathematics, and he is an early contributor to the ADM. Dr. Rach’s fundamental theorems established the basis for the early development of ADM. He was also the first to propose the modified decomposition method (MDM for short). His past and current research continue to advance ADM/MDM work. He prepared the comprehensive bibliography on ADM:
R. Rach, A bibliography of the theory and applications of the Adomian decomposition method, 1961-2011, Kybernetes, Volume 41, Nos. 7/8, 2012, 1087--1148.

Dr. Sergio E. Serrano was born in Santander (Colombia) in 1953. He is a professor of environmental engineering, hydrologic science, applied mathematics, and philosophy at Temple University in Philadelphia. He received his Ph.D. degree at the University of Waterloo (Canada). He has more than one hundred research publications in international science, engineering, and mathematical journals. He is also the author of eight books in environmental engineering, statistics, philosophy, and psychology. Dr. Serrano has been awarded four times with nationally-competitive research grants by the National Science Foundation, Washington, DC. Using adaptations and modifications of the ADM, he has developed hundreds of practical engineering models of flood wave propagation, contaminant transport, and groundwater flow in irregular geometries. He published influential textbooks on applications of differential equations in CAS Maple and Hydrology:
• Sergio E. Serrano, Differential Equations, HydroScience Inc. Amber, Pennsylvania, 2016, ISBN 97809888865211
• Sergio E. Serrano, Hydrology for Engineers, Geologists, and Environmental Professionals. An Integrated Treatment of Surface, Subsurface, and Contaminant Hydrology. HydroScience Inc., Ambler, PA.

Dr. Serrano has a passion for hiking. In 1973, he explored the rain forest of the Darien Sierra on foot from Acandi (Colombia) to Yaviza (Panama). He plays the recorder (Renaissance flute). He enjoys alchemy, archaeology, home wine making, herbal medicine, and cooking. He believes that the joy of meaningful living and meaningful relationships can be found in the simplicity of everyday life. He lives in Philadelphia with his wife of thirty years and his daughter.

Dr. Abdul-Majid Wazwaz is a Professor of Mathematics at Saint Xavier University in Chicago. He received his Ph.D. from the University of Illinois at Chicago. He was the author and co-author of more than four hundred and fifty papers in applied mathematics and mathematical physics. He is the author of five books on the subjects of discrete mathematics, integral equations, and partial differential equations. He has contributed extensively to theoretical advances in solitary waves theory, the ADM, and other computational methods. He is a member of the editorial board of the journals Nonlinear Dynamics (Springer) and Physica Scripta (IOP). For three years in a row, 2014, 2015, and 2016, Thomas Reuters granted him three different badges for being a “Highly Cited Researcher.”

The Adomian decomposition method (ADM for short) has led to several modifications on the method made by various researchers in an attempt to improve the accuracy or expand the application of the original approach. Obviously, this tutorial cannot cover and explain all available improvements for the method. Instead, we show how the Adomian's decomposition method works in a series of examples. The basic spirit of the decomposition method consists of three steps.

• The first one includes a decomposition of the true solution to the initial value problem (or boundary value problem) into the infinite sum of the solution's components to auxiliary problems.
• The second step asks to find the solution's components recursively based on previously obtained components, so it involves solution of a full-history recurrence. This stage requires solving very basic differential equations with explicitly determined right-hand side expressions that are obtained by decomposition of nonlinear terms into the so called Adomian's polynomials.
• Finally, summation of the first finite number of components leads to an approximation to the true solution with any desired level of high precision.

If the given differential equation is homogeneous (without a driving term), then the ADM finite sum always provides the truncated series version to the true solution provided that the Adomian series converges. However, for inhomogeneous differential equations, the ADM finite sum approximation may include some noise terms that eventually are canceled out with subsequent iterations assuming that all iterations are complete. This phenomenon is well-known for Picard's iteration. Unwanted terms in the ADM applications were first discovered by G. Adomian and R. Rach in 1986 (The noisy convergence phenomena in decomposition method solutions, Journal of Computational and Applied Mathematics, Volume 15, Issue 3, July 1986, pages 379--381). Although it is possible to reduce these unwanted noise terms, as shown by Abdul-Majid Wazwaz and other researchers, the issue remains. However, these noise terms usually do not effect the approximations because they are rapidly damped out numerically.

The first proof of convergence of the ADM was given by Cherruault in 1989, who used fixed point theorems for abstract functional equations. Since conditions of these theorems are too restrictive for most physical and engineering applications to be verified in practice, many other articles on convergence of the ADM were published. In spite of the variety of publications on convergence, computational complexity, improvements, and applications of the ADM, no precise criterion of convergence was formulated in the literature, at least in the context of initial-value problems for ordinary and partial differential equations. When the slope function is a holomorphic function, then the formalism of the Cauchy--Kowalevskaya theorem guarantees that solutions of initial-value problems for systems of ordinary differential equations exist and are analytic for small time intervals. Recent reviews on this topic can be found in the following articles: Ray's (New Approach for General Convergence of the Adomian Decomposition Method, World Applied Sciences Journal, 2014, Vol. 32, No. 11, pp. 2264--2268), Rach's (Bougoffa, L., Rach, R.C., El-Manouni, S., A convergence analysis of the Adomian decomposition method for an abstract Cauchy problem of a system of first-order nonlinear differential equations, International Journal of Computer Mathematics, 2013, Volume 90, Issue 2, pages 360--375), and Sunday's (Convergence analysis and implementation of Adomian decomposition method on second-order oscillatory problems, Asian Research Journal of Mathematics, 2017, Vol. 2, No. 5, pp. 1--12. Article no.ARJOM.32011) articles.

Let us compare the ADM and Picard's iteration scheme. The latter can be successfully applied only for differential equations with polynomial slope functions. The ADM extends Picard's requirements for analytically defined nonlinearities. The advantage of Adomian’s method over the Picard scheme is the ease of computation of successive terms. At each iterative step, the Adomian decomposition method actually requires solving the same (very simple) initial value problem with homogeneous initial conditions. Since the ADM presents the solution as (infinite) series, it avoids repetitions that slow down Picard's procedure. Although it is a challenging problem to explicitly determine the radius of convergence for the solution obtained on the basis of the ADM, it is usually much larger than that of Picard's prediction.

These advantages of the ADM over Picard's iteration scheme are abated by preprocessing: one should evaluate Adomian's polynomials that are used at every iteration step. The labor involved in such evaluation grows exponentially in general. This drawback becomes negligible in an explicit one step numerical implementation based on the ADM when only a few first terms are taken into account. It should be noted that the amount of work for the ADM preprocessing is larger than, say, the well known Runge--Kutta or cubic spline algorithms where a user enters only the slope functions and the initial conditions. Usually, the labor required by ADM preprocessing is compensated by a larger step size in the computations than in standard one step algorithms.

We start a demonstration of the Adomian decomposition method with the following initial value problem:

$y' = f(x,y) + g(x), \qquad y(x_0 ) = y_0 ,$
where f is the given (smooth) function, g is an input (driving) term, y is the (unknown) output of the system, and constants x0 and y0 are prescribed. According to the ADM, we seek an approximation to the solution of this initial value problem as an infinite series of functions un to be determined recursively:
$y(x) = \sum_{n\ge 0} u_n = u_0 + u_1 + u_2 + \cdots ,$
while the nonlinear analytic function f(x,y) is assumed to be expanded into an Adomian series of polynomials
$f(x,y) = \sum_{n\ge 0} A_n (x) = A_0 (u_0 ) + A_1 (u_0 , u_1 ) + A_2 \left( u_0 , u_1 , u_2 \right) + \cdots ,$
where An are the (classical) Adomian’s polynomials defined recursively according to the formula:
$A_n \left( u_0 , u_1 , \ldots , u_n \right) = \frac{1}{n!} \, \frac{{\text d}^n}{{\text d} \lambda^n} \left[ f \left( x, \sum_{k\ge 0} u_k \lambda^k \right) \right]_{\lambda =0} , \qquad n=0,1,2,\ldots ;$
which is actually the Faà di Bruno's formula applied to the composition of two functions: f and $$\sum_{k\ge 0} u_k .$$ These polynomials were first introduced by G. Adomian and R. Rach in their 1983 paper (Inversion of nonlinear stochastic operators, Journal of Mathematical Analysis and Applications, Vol. 91, No. 1, January 1983, pp. 39--46). In the preceding formula, the parameter λ is a grouping parameter that is used to identify a certain iteration term. The coefficient of λn corresponds to the nth iteration term. One of the ways to extract this term is to differentiate n times the composition $$f \left( x, \sum_{k\ge 0} u_k \lambda^k \right)$$ and set λ = 0 to eliminate the rest of the terms. Note that the previous terms are eliminated by differentiation. Fortunately, Mathematica has two dedicated commands for extracting coefficients from a polynomial: Coefficient that gives a particular coefficient and CoefficientList that provides all coefficients of the given polynomial (but not an arbitrary function).
The Adomian decomposition method assumes that the solution to the given initial value problem
$\frac{{\text d}y}{{\text d}x} = f(x,y), \qquad y(x_0 ) = y_0 ,$
$y(x) = u_0 + \sum_{n\ge 1} u_n .$
It is also assumed that the slope function can be expanded as an infinite series,
$f(x,y) = \sum_{n\ge 0} A_n \left( u_0 , u_1 , \ldots , u_n \right) ,$
where the An are Adomian's polynomials. To derive the expression for these polynomials, we consider the Taylor expansion of the slope function about u0:
$f(x,y) = f(u_0 ) + f' (u_0 ) \left( y - u_0 \right) + \frac{1}{2!}\, f'' (u_0 ) \left( y - u_0 \right)^2 + \frac{1}{3!}\, f^{(3)} (u_0 ) \left( y - u_0 \right)^3 + \cdots + \frac{1}{n!}\, f^{(n)} (u_0 ) \left( y - u_0 \right)^n + \cdots ,$
where we dropped writing the independent variable for simplicity because it does not affect the derivation. Substituting this expansion into the solution series and noting that
$y - u_0 = \left( u_0 + u_1 + u_2 + \cdots \right) - u_0 = u_1 + u_2 + u_3 + \cdots ,$
yields
\begin{align*} f(y) &= f(u_0 ) + f' (u_0 ) \left( u_1 + u_2 + u_3 + \cdots \right) \\ & \quad + \frac{1}{2!}\, f'' (u_0 ) \left( u_1 + u_2 + u_3 + \cdots \right)^2 + \frac{1}{3!}\, f'' (u_0 ) \left( u_1 + u_2 + u_3 + \cdots \right)^3 + \cdots . \end{align*}
Expanding the first few terms of the Taylor expansion, we get
\begin{align*} \left( u_1 + u_2 + u_3 + \cdots \right)^2 &= u_1^2 + 2u_1 u_2 + 2u_1 u_3 + u_2^2 + 2u_1 u_3 + \cdots , \\ \left( u_1 + u_2 + u_3 + \cdots \right)^3 &= u_1^3 + 3u_1^2 u_2 + 3u_1^2 u_3 + 3u_1 u_2^2 + 6u_1 u_2 u_3 + u_2^3 + \cdots . \end{align*}
Substituting these expansions into the Taylor expansion of the slope function gives
\begin{align*} f(y) &= f(u_0 ) + f' (u_0 ) \left( u_1 + u_2 + u_3 + \cdots \right) \\ & \quad + \frac{1}{2!}\, f'' (u_0 ) \left( u_1^2 + 2u_1 u_2 + 2 u_1 u_3 + u_2^2 + 2 u_1 u_4 + 2 u_2 u_3 + \cdots \right) \\ & \quad + \frac{1}{3!}\, f''' (u_0 ) \left( u_1^3 + 3u_1^2 u_2 + 3 u_1^2 u_3 + 3 u_1 u_2^2 + 6u_1 u_2 u_3 + u_2^3 + \cdots \right) + \cdots . \end{align*}
After rearranging the terms,
\begin{align*} f(y) &= f(u_0 ) + f' (u_0 )\, u_1 + f' (u_0 ) \,u_2 + \frac{1}{2!}\, f'' (u_0 ) \, u_1^2 + f' (u_0 )\, u_3 + \frac{1}{2!}\, f'' (u_0 ) \, 2u_1 u_2 \\ & \quad + \frac{1}{3!}\, f''' (u_0 ) \, u_1^3 + f' (u_0 )\, u_4 + \frac{1}{2!}\, f'' (u_0 ) \left( 2 u_1 u_3 + u_2^2 \right) + \frac{1}{4!}\, f^{(4)} (u_0 ) \, u_1^4 + \cdots . \end{align*}
To calculate the first Adomian polynomial A0, consider any term where the sum of the subscripts of the ui add up 0. Next to get A1, consider each term where the subscripts of the ui add up 1. In a similar fashion, the general Adomian polynomial An includes any term where the sum of the subscripts of the ui add up to n.It is worth noting that the sum of the subscripts of the term uij is ij and not i.    ▣

There are several known modifications of the Adomian polynomials, but we will use only one, known as accelerated Adomian polynomials (since there is no standard notation, we denote them by aAn) that were first introduced by George Adomian in his 1989 book (Nonlinear Stochastic Systems Theory and Applications to Physics, Kluwer Academic Publishers, Dordrecht):
$aA_n \left( u_0 , u_1 , \ldots , u_n \right) = f\left[ u_0 + u_1 + \cdots + u_n \right] - f\left[ u_0 + u_1 + \cdots + u_{n-1} \right] , \qquad J_0 = f\left[ u_0 \right] .$

In 2000, Khelifa and Cherruault established a bound for Adomian polynomials (Kybernetes, Vol. 29, No. 3, pp. 332--355):

$\left\vert A_n \right\vert \le \frac{(n+1)^n}{(n+1)!} \, M^{n+1} , \qquad \sup \left\vert f^{(k)} (x,y_0 ) \right\vert \le M .$
In general, since the Adomian method requires the existence of all derivatives (holomorphy) for the slope function, which is more restrictive than the Lipschitz condition required for Picard's procedure, it is expected that the ADM converges faster than its Picard counterpart.

The initial term u0 in the series expansion $$y(x) = \sum_{n \ge 0} u_n (x)$$ is a solution of the truncated initial value problem:

$y' = g(x), \qquad y(x_0 ) = y_0 ,$
because it incorporates the influence of all inhomogeneous terms of the problem (that includes the initial condition and the forcing input function g). Once it is achieved, all other terms in the infinite series are solutions of the homogeneous initial value problems. Therefore, u0 is usually not the same initial approximation as in Picard's iteration procedure that is just the initial value y0. It should be noted that when the Adomian method is applied to other problems (not necessarily initial value problems for first order differential equations as we consider in this section) one can choose another approximation to the nonhomgeneous differential equation: it is important to exclude, depending on x, the input term from the given slope function. So in this section, we follow the standard ADM procedure applicable to differential equations of the form: $$y' = f(x,y) + g(x) .$$ Therefore, the initial term u0 could be chosen with some flexibility. Once u0 is determined, all other terms in Adomian's decomposition follow from the recursive initial value problems of Appell's type:
$\frac{{\text d}u_n}{{\text d}x} = A_{n-1} \left( u_0 , u_1 , \ldots , u_{n-1} \right) , \qquad u_n (x_0 ) = 0 , \qquad n=1,2,\ldots ,$
where the Adomian polynomials are determined recursively
\begin{eqnarray*} A_0 &=& f(x, u_0 ) , \\ A_1 &=& u_1 f' (u_0 ) , \\ A_2 &=& u_2 f' (u_0 ) + \frac{1}{2!}\, u_1^2 f'' (u_0 ) , \\ A_3 &=& u_3 f' (u_0 ) + \frac{1}{3!}\, u_1^3 f''' (u_0 ) + \frac{2}{2!} \,u_1 u_2 f'' (u_0 ) , \\ A_4 &=& u_4 f' (u_0 ) + \left[ \frac{1}{2!}\, u_2^2 + u_1 u_3 \right] f'' (u_0 ) + \frac{1}{2!} \,u_1^2 u_2 f''' (u_0 ) + \frac{1}{4!} \, u_1^4 f'''' (u_0 ) , \\ A_5 &=& u_5 f' (u_0 ) + f'' (u_0 ) \left( u_1 u_4 + u_2 u_3 \right) + \frac{1}{2!}\,f''' (u_0 ) \left( u_1^2 u_3 + u_1 u_2^2 \right) + \frac{1}{3!}\, f^{(4)} (u_0 ) \, u_1^3 u_2 + \frac{1}{5!} \, f^{(5)} (u_0 ) \, u_1^5 . \end{eqnarray*}
Determination of other Adomian polynomials requires more work.
\begin{eqnarray*} A_6 &=& u_6 f' (u_0 ) + f'' (u_0 ) \left( u_1 u_5 + u_2 u_4 + u_3^2 \right) + \frac{1}{3!}\,f''' (u_0 ) \left( 3\,u_1^2 u_4 + u_2^3 + 6\,u_1 u_2 u_3 \right) \\ &\quad& + \frac{1}{4!}\,f^{(4)} (u_0 ) \left( 4\,u_1^3 + 6\,u_1^2 u_2^2 \right) + \frac{1}{4!}\,f^{(5)} (u_0 )\, u_1^4 u_2 + \frac{1}{6!} \,f^{(6)} (u_0 ) \, u_1^6 , \\ A_7 &=& u_7 f' (u_0 ) + f'' (u_0 ) \left( u_1 u_6 + u_2 u_5 + u_3 u_4 \right) + \frac{1}{2}\,f''' (u_0 ) \left( u_1^2 u_5 + u_1 u_3^2 + u_3 u_2^2 + 2\,u_1 u_2 u_4 \right) \\ & \quad& + \frac{1}{3!} \,f^{(4)} (u_0 ) \left( u_1^3 u_4 +3\, u_1^2 u_2 u_3 + u_1 u_2^3 \right) + \frac{1}{4!}\,f^{(5)} (u_0 ) \left( u_1^4 u_3 + 2\,u_1^3 u_2^2 \right) + \frac{1}{6!} \,f^{(6)} (u_0 ) \,u_1^5 u_2 + \frac{1}{7!} \,f^{(7)} (u_0 ) \,u_1^7 , \\ A_8 &=& u_8 f' (u_0 ) + f'' (u_0 ) \left( u_1 u_7 + u_2 u_6 + u_3 u_5 + u_4^2 /2 \right) \\ & \quad& + \frac{1}{2}\,f''' (u_0 ) \left( u_1^2 u_6 + u_2^2 u_4 + u_2 u_3^2 + 2\,u_1 u_2 u_5 + 2\, u_1 u_3 u_4 \right) + \frac{1}{3!} \,f^{(4)} (u_0 ) \left( u_1^3 u_5 +3\, u_1^2 u_2 u_4 + 3\, u_1 u_2^2 u_3 + 3\,u_1^2 u_3^2 /2 + u_2^4 /4 \right) \\ &\quad& + \frac{1}{4!}\,f^{(5)} (u_0 ) \left( u_1^4 u_4 + 4\,u_1^3 u_2 u_3 + 2\, u_1^2 u_2^3 \right) + \frac{1}{6!} \,f^{(6)} (u_0 ) \left( u_1^5 u_3 + 15\, u_1^4 u_2^2 \right) + \frac{1}{6!} \,f^{(7)} (u_0 ) \,u_1^6 u_2 + \frac{1}{8!} \,f^{(8)} (u_0 ) \, u_1^8 , \\ \vdots &=& \vdots . \end{eqnarray*}
Note: Adomian polynomials are not polynomials in the variable x. They are only polynomials in u1, … ,un, but not in u0. Additionally, components u1, … ,un might not be polynomials in the variable x.

Adomian's polynomials can be evaluated using the Duan's code (Applied Mathematical Modelling, 2013, 37, Issues 21-22, pp. 8687--8708) as follows:

A[0]=f[u[0]];
For[n = 1,n<=M,n++,c[n,1]=u[n];
For[k = 2,k<=n,k++,
c[n,k]=Expand[1/n*Sum[(j + 1)*u[j + 1]* c[n-1-j,k-1],{j,0,n-k}]]];
der = Table[D[f[u[0]],{u[0],k}],{k,1,M}];
A[n]=Take[der,n].Table[c[n,k],{k,1,n}]];
Table[A[n],{n,0,M}]]

Example: Let us consider the nonlinear (pendulum) operator: $$N\left[ \theta \right] = \sin\theta .$$ Then the first few Adomian's polynomials for this nonlinear operator are

\begin{eqnarray*} A_0 &=& \sin \theta_0 , \\ A_1 &=& \theta_1 \cos \theta_0 , \\ A_2 &=& - \frac{1}{2} \,\theta_1^2 \sin \theta_0 + \theta_2 \cos \theta_0 , \\ A_3 &=& - \frac{1}{6}\, \theta_1^3 \cos \theta_0 - \theta_1 \theta_2 \sin \theta_0 + \theta_3 \cos \theta_0 , \end{eqnarray*}
and so on. The modified polynomials are
\begin{eqnarray*} J_1 &=& \sin \left( u_0 + u_1 \right) - \sin u_0 , \\ J_2 &=& \sin \left( u_0 + u_1 + u_2 \right) - \sin \left( u_0 + u_1 \right) , \end{eqnarray*}
and so on. Let us consider the hyperbolic sine nonlinear function $$N\left[ x \right] = \sinh \left( x/2 \right) .$$ The corresponding Adomian's polynomials are
\begin{eqnarray*} A_1 &=& \frac{1}{2}\, x_1 \cosh \left( x_0 /2 \right) , \\ A_2 &=& \frac{1}{2}\, x_2 \cosh \left( \frac{x_0}{2} \right) + \frac{1}{2} \, x_1^2 \frac{1}{2^2} \,\sinh \left( \frac{x_0}{2} \right) , \\ A_3 &=& \frac{1}{3}\, x_3 \cosh \left( \frac{x_0}{2} \right) + \frac{1}{2^2} x_1 x_2 \sinh \left( \frac{x_0}{2} \right) + \frac{1}{6} \,x_1^3 \frac{1}{2^3} \,\cosh \left( \frac{x_0}{2} \right) \\ A_4 &=& \frac{1}{2}\, x_4 \cosh \left( \frac{x_0}{2} \right) + \left[ \frac{1}{2}\, x_2^2 + x_1 x_2 \right] \frac{1}{2^2} \sinh \left( \frac{x_0}{2} \right) + \frac{1}{2}\, x_1^2 x_2 \frac{1}{2^3} \,\cosh \left( \frac{x_0}{2} \right) + \frac{1}{24}\, \frac{1}{2^4} \,x_1^4 \sinh \left( \frac{x_0}{2} \right) , \end{eqnarray*}
and so on. Our next example is the power function: $$N\left[ u \right] = u^{\gamma} ,$$ where γ is a real number. Then corresponding Adomian's polynomials are
\begin{eqnarray*} A_0 &=& u_0^{\gamma} , \\ A_1 &=& \gamma \,u_0^{\gamma -1} u_1 , \\ A_2 &=& \gamma \,u_0^{\gamma -1} u_2 + \frac{1}{2}\,\gamma \left( \gamma -1 \right) u_0^{\gamma -2} u_1^2 , \\ A_3 &=& \gamma \,u_0^{\gamma -1} u_3 + \gamma \left( \gamma -1 \right) u_0^{\gamma -2} u_1 u_2 + \frac{1}{6}\, \gamma \left( \gamma -1 \right) \left( \gamma -2 \right) u_0^{\gamma -3} u_1^3 , \\ A_4 &=& \gamma \,u_0^{\gamma -1} u_4 + \gamma \left( \gamma -1 \right) u_0^{\gamma -2} \left[ \frac{1}{2}\, u_2^2 + u_1 u_3 \right] + \frac{1}{2} \, \gamma \left( \gamma -1 \right) \left( \gamma -2 \right) u_0^{\gamma -3} u_1^2 u_2 + \frac{1}{24} \, \gamma \left( \gamma -1 \right) \left( \gamma -2 \right) \left( \gamma -3 \right) u_0^{\gamma -4} u_1^4 . \end{eqnarray*}
■

Example: Consider the simplest algebraic equation---the quadratic equation

$a\,x^2 + b\,x + c = 0, \qquad b\ne 0,$
where 𝑎, b, and c are some given constants. The case b = 0 can be easily included in our consideration by adding and subtracting the x term. The idea of isolating x is to transform the given problem to the form x = Φ(x), with some continuous function Φ, so that it appears as a fixed point problem. From the given quadratic equation, we obtain
$x = - \frac{c}{b} - \frac{a}{b}\, x^2 = \alpha + \gamma\,x^2 \qquad\mbox{with } \quad \alpha = - \frac{c}{b}, \ \gamma = - \frac{a}{b} .$
In general, one can add and subtract k x, k ≠ 0, and express x as
$x = -\frac{c}{k} + \frac{k-b}{k}\, x - \frac{a}{k}\, x^2 = \alpha + \beta\, x + \gamma \,x^2 \qquad (k \ne 0).$
In our case, we have k = b, but later we will consider other values of k. It is helpful to work a numerical example and show all calculations for it. Hence, we consider the equation
$x^2 - x -6 =0,$
whose solutions are known to be x = 3 and x = -2. Let us see how ADM arrives at these values. We rewrite our quadratic equation in the form suitable for the fixed point theorem
$x = -6 + x^2 \qquad\Longrightarrow \qquad \alpha = -6, \ \beta = 0, \ \gamma = 1.$
We decompose its solution into an infinite series
$x = x_0 + x_1 + x_2 + \cdots = \sum_{n\ge 0} x_n .$
The quadratic equation can rewritten in the operator form
$L \,x = N\left[ x \right] + g, \qquad L\,x = x, \quad N\left[ x \right] = \beta\,x + \gamma\,x^2 = x^2 , \quad g = \alpha = -6.$
Our next step is to decompose the nonlinear term x² as a series over Adomian's polynomials:
$N\left[ x \right] = \beta\,x + \gamma\,x^2 = \sum_{n\ge 0} A_n .$
The components xn of the decomposition solution $$x= \sum_{n\ge 0} x_n$$ are determined recursively
\begin{align*} x_0 &= - \frac{c}{b} = -\frac{c}{k} = \alpha \qquad\Longrightarrow \qquad x_0 = -6, \\ x_{n+1} &= A_n , \qquad n=0,1,2,\ldots . \end{align*}
$\frac{-\beta - \beta \sqrt{1 - 4\alpha \gamma /\beta^2}}{2\gamma} =$
The series converges when $$\left\vert 4\alpha \gamma /\beta^2 \right\vert < 1 .$$ So the series converges when $$b^2 > |a\,c| .$$ This condition is violated for our case and the Adomian series diverges. We can see this by finding derivatives of f(x) = x² - 5 at the roots: $$f' (3) = 2\,(3) = 6 > 1$$ and $$f' (-2) = -4 ,$$ whose absolute value exceeds 1. Therefore, our fixed point equation x = f(x) is not a contractor (as required by the condition discovered by Cherruault for the ADM convergence).

To overcome this misfortune, we add and subtract 4x to the both sides of the given equation

$x^2 -x -6 + 4x = 4x \qquad \Longrightarrow \qquad x = f(x) = \frac{1}{4}\, x^2 + \frac{3}{4}\,x - \frac{3}{2} = N[x] - \frac{3}{2} .$
Now we choose the initial approximation to be x0 = -3/2, and expand the remaining function N[x] = x²/4 - 3x/4 into the Adomian series
$\frac{1}{4}\, x^2 + \frac{3}{4}\,x = \frac{x}{4} \left( x+3 \right) = \sum_{n\ge 0} A_n ,$
where N[x] = x(x + 3)/4. Substituting the series into the quadratic equation, we obtain
$x_{n+1} = A_n , \qquad n=0,1,2,\ldots .$
Now we calculate the Adomian polynomials using derivatives
$N(x) = \frac{x}{4} \left( x+3 \right) , \qquad N'(x) = \frac{2x+3}{4}, \qquad N''(x) = \frac{1}{2} .$
Therefore, \begin{align*} x_1 &= A_0 = N\left( -\frac{3}{2} \right) = - \frac{9}{16} = -0.5625, \\ x_2 &= A_1 = x_1 N'\left( -\frac{3}{2} \right) = 0, \\ x_3 &= A_2 = x_2 N'\left( -\frac{3}{2} \right) + x_1^2 N'' /2 = \frac{81}{1024} = 0.0791016 , \\ x_4 &= A_3 = x_3 N' \left( -\frac{3}{2} \right) + x_1 x_2 /2 = 0, \\ x_5 &= A_4 = (x_2^2 + 2\,x_1 x_3)/4 = - \frac{729}{32 768} \approx -0.0222473 , \\ x_6 &= A_5 = 0, \\ x_7 &= A_6 = \frac{19 683}{2 097 152} \approx 0.00938559 . \end{align*}
f[x_] = x*(x + 3)/4
f1[x_] = (2*x + 3)/4
x0 = -3/2
x1 = f[x0]
x2 = x1*f1[x0]
x3 = x2*f1[x0] + x1^2/4
x4=x3*f1[x0]+x1*x2/2
x5=(x2^2 + 2*x1*x3)/4
x6=0
x7=(x1 x5 + x2 x4 + x3^2)/2
xx7=x0+x1+x3+x5+x7
-(4186461/2097152)
Adding the first seven terms, we obtain the approximation of the root to be -1.99626. However, if we add one more term, we get a better approximation -1.99978 to the actual root. Mathematica confirms
x9 = (x1 x7 + x3 x5)/2
x9 + xx7
-(33550737/16777216)

Now we find another root by choosing another modification:
$x^2 -x -6 - 4x = -4x \qquad \Longrightarrow \qquad x = f(x) = -\frac{1}{4}\, x^2 + \frac{5}{4}\,x + \frac{3}{2} = N[x] + \frac{3}{2} .$
Using N[x] = -x(x - 5)/4, we calculate the Adomian polynomials: \begin{align*} x_0 &= \frac{3}{2} , \\ x_1 &= A_0 = N\left( \frac{3}{2} \right) = \frac{21}{16} = 1.3125 , \\ x_2 &= A_1 = x_1 N'\left( \frac{3}{2} \right) = \frac{21}{32} \approx 0.65625 , \\ x_3 &= A_2 = x_2 N'\left( \frac{3}{2} \right) + x_1^2 N'' /2 = - \frac{105}{1024} \approx -0.102539 , \\ x_4 &= A_3 = x_3 N'\left( \frac{3}{2} \right) -x_1 x_2 /2 = - \frac{987}{2038} \approx -0.481934 . \end{align*} When we add the first seven components, we get 2.87075.
f[x_] = -x*(x - 5)/4
f1[x_] = (-2*x + 5)/4
x0 = 3/2
x1 = f[x0]
x2 = x1*f1[x0]
x3 = x2*f1[x0] - x1^2/4
x4 = x3*f1[x0] - x1*x2/2
x5 = x4*f1[x0] - (x2^2 + 2*x1*x3)/4
x6 = x5*f1[x0] - (x1*x4 + x2*x3)/2
x7 = x6*f1[x0] - (x1 x5 + x2 x4 + x3^2)/2
xx7 = x0 + x1 + x3 + x5 + x7
6020397/2097152
■

Example: Consider the initial value problem for the linear equation

$y' = x \left( 3 - 2y \right) = f(x,y) + g(x), \qquad y(0) =1,$
where $$g(x) = 3x$$ and $$f(x,y) = -2x\,y .$$ In this example, we will find the first 5 terms of the Adomian approximation for this IVP. We will then compare our approximation with the true solution found with Mathematica
DSolve[{y'[x] == x*(3 - 2*y[x]), y[0] == 1}, y[x], x]// Flatten
{y[x] -> 1/2 E^-x^2 (-1 + 3 E^x^2)}
$y= \phi (x) = \frac{3}{2} - \frac{1}{2}\, e^{-x^2} = 1 + \frac{x^2}{2} - \frac{x^4}{4} + \frac{x^6}{12} - \frac{x^8}{2\cdot 4!} + \frac{x^{10}}{2\cdot 5!} + \cdots ,$
in order to see how good our approximation is. According to Adomian, we seek its solution as an infinite sum (which we truncate to keep the first few terms)
$y(x) = \sum_{n\ge 0} u_n (x) = u_0 (x) + u_1 (x) + u_2 (x) + \cdots ,$
where the first term u0 is the solution to the following initial value problem:
$u' = g(x) = 3x , \qquad u(0) =1 ,$
because we need to start with the nonhomogeneous part of the original initial value problem. Since this is a calculus problem, its solution is $$u_0 (x) = C + \frac{3\,x^2}{2} ,$$ where C is a constant of integration which happens to be the initial condition of the IVP. Hence, C=1. This u0 differs from Picard's starting term, which is 1. With u0, we can find A0, the zeroth Adomian polynomial:
$A_0 = f(x, u_0 ) = -2x \,u_0 = -2x - 3\,x^3 .$
f[x_, y_] = -2*x*y
f1[x_] = -2*x
u0[x_] = 1 + 3*x^2/2
A0[x_] = f[x, u0[x]]
-2 x (1 + (3 x^2)/2)
All Adomian polynomials, except the first one, are computed according to a single formula (which is valid only for linear equations)
$A_n (x) = -2x\,u_n , \qquad n=1,2,\ldots .$
Mathematica confirms:
f[u_] := - 2*u*x
{-2 x u[0], -2 x u[1], -2 x u[2], -2 x u[3], -2 x u[4], -2 x u[5]}
Also, we can achieve this with another Mathematica command (which can only be applied for polynomial input):
CoefficientList[-2 x*(u0 + u1*s + u2*s^2 + u3*s^3 + u4*s^4 + u5*s^5), s, 5]
{ -2 u0 x, -2 u1 x, -2 u2 x, -2 u3 x, -2 u4 x}
However, we don't need Adomian polynomials in our simple linear differential equation. We just substitute the infinite series form $$y(x) = \sum_{n\ge 0} u_n (x)$$ into the given equation:
$\frac{\text d}{{\text d}x} \left[ u_0 + u_1 + u_2 + \cdots \right] = 3x - 2x \left[ u_0 + u_1 + u_2 + \cdots \right] .$
This leads to a sequence of initial value problems:
$\frac{\text d}{{\text d}x} \, u_{n+1} = - 2x \, u_n , \quad u_{n+1} (0) = 0, \qquad n=0,1,2,\ldots .$
This differential equation is similar to the recursive equation for Appell polynomials. We ask Mathematica to find the solution to the above IVP.
DSolve[{y'[x] == -2*x*(1 +3* x^2 /2), y[0] == 0}, y[x], x]// Flatten
u[1][x_] = y[x] /. %
{y[x] -> 1/4 (-4 x^2 - 3 x^4)}
Hence, $$u_1 = - x^2 - 3\,x^4 /4$$ and we get the first order approximation
$y_1 = u_0 + u_1 = 1 + \frac{3\,x^2}{2} - x^2 - \frac{3\,x^4}{4} = 1 + \frac{1}{2}\,x^2 - \frac{3\,x^4}{4} ,$
phi1[x_] = Simplify[1 + 3*x^2/2 + u[1][x]]
1 + x^2 /2 - 3 x^4 /4
which differs from Picard's, $$\phi_2 = 1 + \frac{x^2}{2} - \frac{x^4}{4} ,$$ in the coefficient of the largest power of x. Next, we find u2 by solving the initial value problem
$y' = A_1 = -2x \, u_1 = 2\,x^3 + 3\,x^5 /2, \quad y(0) =0 .$
So we ask Mathematica again to do this job for us:
DSolve[{y'[x] == -2*x*u[1][x], y[0] == 0}, y[x], x] //Flatten
{y[x] -> 1/4 (2 x^4 + x^6)}
u[2][x_] = y[x] /. %
Adding to the previous sum, we get the next approximation
phi2[x_] = Simplify[phi1[x] + u[2][x]]
1/4 (4 + 2 x^2 - x^4 + x^6)
Similarly,
DSolve[{y'[x] == -2*x*u[2][x], y[0] == 0}, y[x], x] //Flatten
u[3][x_] = y[x] /. %
{y[x] -> 1/48 (-8 x^6 - 3 x^8)}
So we get
phi3[x_] = Simplify[phi2[x] + u[2][x]]
1/48 (48 + 24 x^2 - 12 x^4 + 4 x^6 - 3 x^8)
$u_0 = 1 + \frac{3\,x^2}{2} , \quad u_1 = -x^2 - \frac{3}{4}\, x^4 , \quad u_2 = \frac{1}{2} \,x^4 + \frac{1}{4}\, x^6 , \quad u_3 = -\frac{1}{6}\,x^6 - \frac{1}{16}\, x^8 , \quad u_4 = \frac{1}{24} \,x^8 + \frac{1}{80}\,x^{10} , \quad u_5 = -\frac{1}{120} \,x^{10} - \frac{1}{480}\, x^{12} .$
$\phi_5 (x) = u_0 + u_1 + u_2 +u_3 +u_4 +u_5 = 1+ \frac{x^2}{2} - \frac{x^4}{4} + \frac{1}{12}\,x^6 - \frac{1}{48}\, x^8 + \frac{1}{240}\, x^{10} - \frac{1}{480}\, x^{12} .$
Now we plot the exact solution along with our 5-th order Adomian approximation.
Plot[{Callout[3/2-Exp[-x^2]/2, "exact"], Callout[phi5[x], "ADM approximation"]}, {x,0,2.0}, PlotStyle->Thick]
■

Example: Consider the initial value problem for the Riccati equation

$y' = x^2 + y^2 , \qquad y(0) =1.$
According to Mathematica, its explicit solution is known to be (see Riccati section)
riccati = DSolve[{y'[x] == x^2 + (y[x])^2, y[0] == 1}, y[x], x]; phi[x_] = y[x] /. riccati
$y(x) = \frac{x^2 J_{5/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{3}{4} \right) - x^2 J_{1/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{1}{4} \right) + J_{1/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{3}{4} \right) - x^2 J_{3/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{3}{4} \right) }{x\, J_{1/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{1}{4} \right) -2\,J_{-1/4} \left( \frac{x^2}{2} \right) \Gamma \left( \frac{3}{4} \right) } ,$
DSolve[{y'[x] == x^2 + (y[x])^2, y[0] == 1}, y[x], x]
where Jν(z) is the Bessel function of the first kind and $$\Gamma (z) = \int_0^{\infty} e^{-t} t^{z-1}\,{\text d}t$$ is the gamma function of Euler.

This function blows up when the denominator is zero. This point can be determined approximately from the following graph.

y[x_] = (-x^2 BesselJ[-(3/4), x^2/2] Gamma[1/4] + x^2 BesselJ[-(5/4), x^2/2] Gamma[3/4] + BesselJ[-(1/4), x^2/2] Gamma[3/4] - x^2 BesselJ[3/4, x^2/2] Gamma[3/ 4])/(x (BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4]))
deno[x_] = (x (BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4]))
Plot[{y[x], deno[x]}, {x, 0, 3}, PlotStyle -> {{Thick, Blue}, {Thick, Red}}]

By plotting the denominator of Ricatti solution, we find approximately the point where its graph (in red) crosses the horizontal axis. This point indicates the position where the solution to the Riccati equation blows up. It is impossible to determine this point from the initial value problem only. Fortunately, Mathematica can help with its determination.
FindRoot[deno[x], {x, 1}]
{x -> 0.969811}
Therefore, we know now that the solution (which is unique) to the given initial value problem exists only within the interval [0, 0.9698).

Our next step is to determine a power series representation for the solution. We can find its 10-term Maclaurin series expansion with a one line Mathematica command:

AsymptoticDSolveValue[{y'[x] == x^2 + (y[x])^2, y[0] == 1}, y[x], {x, 0, 10}]
1 + x + x^2 + (4 x^3)/3 + (7 x^4)/6 + (6 x^5)/5 + (37 x^6)/30 + ( 404 x^7)/315 + (369 x^8)/280 + (428 x^9)/315 + (1961 x^10)/1400
$y(x) = 1 + x + x^2 + \frac {4} {3}\, x^3 + \frac {7} {6}\, x^4 + \frac{6}{5}\, x^5 + \frac{37}{30}\, x^6 + \frac{404}{315}\, x^7 + \cdots .$
At this point, although we know many terms in its Maclaurin series, it is impossible to find its radius of convergence because we don't know how coefficients grow.

According to G. Adomian, we seek the solution as an infinite sum

$y(x) = \sum_{n\ge 0} u_n ,$
where u0 is the unique solution of the following IVP:
$y' = x^2 , \qquad y(0) =1;$
which includes only the nonhomogeneous components of the original IVP. Using Mathematica, we obtain
DSolve[{y'[x] == x^2 , y[0] == 1}, y[x], x]
u[0][x_] = y[x] /. %
u[0][x_]:= 1/3 (3 + x^3)
So
$u_0 = 1 + \frac{1}{3}\, x^3 .$
Since the problem is nonlinear, we need Adomian's polynomials to expand the nonlinear term (y²):
f[u_] := u*u
{u[0]^2, 2 u[0] u[1], u[1]^2 + 2 u[0] u[2], 2 u[1] u[2] + 2 u[0] u[3], 2 (u[2]^2/2 + u[1] u[3]) + 2 u[0] u[4], 2 (u[2] u[3] + u[1] u[4]) + 2 u[0] u[5]}
We can calculate the Adomian polynomials manually upon substituting the series $$y= \sum_n u_n$$ into the differential equation:
$\frac{{\text d}}{{\text d}x} \left[ u_0 + u_1 + u_2 + \cdots \right] = x^2 + \left[ u_0 + u_1 + u_2 + \cdots \right]^2 = x^2 + A_0 + A_1 + A_2 + \cdots .$
In our case of the standard Riccati equation, we have
\begin{eqnarray*} A_0 &=& u_0^2 = \frac{1}{9} \left( 3 + x^3 \right)^2 , \\ A_1 &=& 2\,u_1 \,u_0 = \frac{2}{3} \left( 3 + x^3 \right) u_1 , \\ A_2 &=& 2\,u_2 \,u_0 + u_1^2 = \frac{2}{3} \left( 3 + x^3 \right) u_2 + u_1^2 , \\ A_3 &=& 2\,u_3 \,u_0 + 2\,u_1 u_2 = \frac{2}{3} \left( 3 + x^3 \right) u_3 + 2\,u_1 u_2 , \\ A_4 &=& 2\,u_4 \,u_0 + u_2^2 + 2\, u_1 u_3 = \frac{2}{3} \left( 3 + x^3 \right) u_4 + u_2^2 + 2\, u_1 u_3 , \\ A_5 &=& 2\,u_5 \,u_0 + 2 \left( u_1 u_4 + u_2 u_3 \right) , \end{eqnarray*}
and so on. Then we calculate components in the Adomian series solution $$y= \sum_n u_n$$ by solving a sequence of initial value problems:
$\frac{\text d}{{\text d}x}\, u_{n+1} = A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (0) =0 .$
Since
$A_0 (x) = \left( 1 + \frac{1}{3}\, x^3 \right)^2 ,$
we find u1 by solving the initial value problem
$u'_1 = A_0 (x) = \left( 1 + \frac{1}{3}\, x^3 \right)^2 , \qquad u_1 (0) =0 .$
A[0][x_] = Expand[(1 + x^3 /3!)^2]
1 + (2 x^3)/3 + x^6/9
DSolve[{y'[x] == A[0][x], y[0] == 0}, y[x], x] // Flatten
u[1][x_] = y[x] /. %
1/126 (126 x + 21 x^4 + 2 x^7)
Next iteration:
A[1][x_] := Simplify[2*u[0][x]*u[1][x]]
2 x + x^4 + x^7/7 + (2 x^10)/189
DSolve[{y'[x] == A[1][x], y[0] == 0}, y[x], x] // Flatten
u[2][x_] := y[x] /. %
(83160 x^2 + 16632 x^5 + 1485 x^8 + 80 x^11)/83160
phi2[x_] := Simplify[u[0][x] + u[1][x] + u[2][x]]
1 + x + x^2 + x^3/3 + x^4/6 + x^5/5 + x^7/63 + x^8/56 + (2 x^11)/2079
Comparing the two-term Adomian approximation φ2(x) with the true Maclaurin expansion, we see that φ2(x) has correct first three terms and the noise terms:
$\phi_2 (x) = 1 + x + x^2 + \underbrace{\frac{x^3}{3} + \frac{x^4}{6} + \frac{x^5}{5} + \frac{x^7}{63} + \frac{x^8}{56} + \frac{2\,x^{11}}{2079}}_{\mbox{noise terms}} .$
The next Adomian polynomial is found with the aid of Mathematica
A2[x_] = Simplify[2*u0[x]*u2[x] + (u1[x])^2]
3 x^2 + (7 x^5)/5 + (8 x^8)/35 + (53 x^11)/2772 + (13 x^14)/14553
DSolve[{y'[x] == A2[x], y[0] == 0}, y[x], x] // Flatten
(3492720 x^3 + 814968 x^6 + 88704 x^9 + 5565 x^12 + 208 x^15)/3492720
u[3][x_] = (3492720 x^3 + 814968 x^6 + 88704 x^9 + 5565 x^12 + 208 x^15)/3492720
phi3[x_] = Simplify[phi2[x] + u3[x]]
1 + x + x^2 + (4 x^3)/3 + x^4/6 + x^5/5 + ( 7 x^6)/30 + x^7/63 + x^8/56 + (8 x^9)/315 + (2 x^11)/2079 + ( 53 x^12)/33264 + (13 x^15)/218295
The three term approximation $$\phi_3 (x) = u_0 + u_1 + u_2 + u_3$$ is
\begin{align*} \phi_3 (x) &= u_0 + u_1 + u_2 + u_3 = 1 + x + x^2 + \frac{4}{3}\, x^3 \\ &\quad + \underbrace{\frac{x^4}{6} + \frac{x^5}{5} + \frac{7}{30}\, x^6 + \frac{x^7}{63} + \frac{1}{56}\, x^8 + \frac{8}{315}\, x^9 + \frac{2}{2079}\, x^{11} + \frac{53}{33264}\, x^{12} + \frac{13}{218295}\,x^{15}}_{\mbox{noise terms}} . \end{align*}
For the next term in Adomian's series, we need to find A3:
A3[x_] = Simplify[2*u0[x]*u3[x] + 2*u1[x]*u2[x]]
4 x^3 + (28 x^6)/15 + (143 x^9)/420 + (17 x^12)/495 + ( 3613 x^15)/1746360 + (46 x^18)/654885
DSolve[{y'[x] == A3[x], y[0] == 0}, y[x], x] // Flatten
u4[x_] = 1/ 103524220800 (103524220800 x^4 + 27606458880 x^7 + 3524753232 x^10 + 273490560 x^13 + 13386165 x^16 + 382720 x^19)
Then we calculate the corresponding sum of the first four terms:
phi4[x_] = Simplify[phi3[x] + u4[x]]
1 + x + x^2 + (4 x^3)/3 + (7 x^4)/6 + x^5/5 + (7 x^6)/30 + ( 89 x^7)/315 + x^8/56 + (8 x^9)/315 + (143 x^10)/4200 + ( 2 x^11)/2079 + (53 x^12)/33264 + (17 x^13)/6435 + ( 13 x^15)/218295 + (3613 x^16)/27941760 + (46 x^19)/12442815
\begin{align*} \phi_4 (x) &= \phi_3 + u_4 = 1 + x + x^2 + \frac{4}{3}\, x^3 + \frac{7}{6}\, x^4 \\ & \quad + \underbrace{\frac{x^5}{5} + \frac{7}{30}\, x^6 + \frac{89}{315}\, x^7 + \frac{x^8}{56} + \frac{8}{315}\, x^9 + \frac{143}{4200}\, x^{10} + \frac{2}{2079} \, x^{11} + \frac{53}{33264} \, x^{12} + \frac{17}{6435}\, x^{13} + \frac{13}{218295}\, x^{15} + \frac{3613}{27941760} \, x^{16} + \frac{46}{12442815} \,x^{19}}_{\mbox{noise terms}} . \end{align*}
This series reveals that the first terms up to the power 4 are correct. We next compare the last Adomian approximation with the true solution by calculating the error function.
Plot[{Callout[y[x], "exact solution"], Callout[phi4[x], "ADM 4-approximation"]}, {x, 0, 0.7}, PlotStyle -> Thick]
error[x_] = y[x] - phi4[x]
Plot[{error[x]}, {x, 0, 0.7}, PlotStyle -> {Thick}]
We compare our ADM approximation with true solution and approximation obtained with the aid of classical Runge--Kutta of order 4 (abbreviated by RK4) with step size h = 0.1.
Point True value ADM-4 approximation RK4 approximation
x = 0.1 1.11146 1.11145 1.11146
x = 0.2 1.25302 1.25262 1.25302
x = 0.3 1.43967 1.43617 1.43967
Expantion of the true solution into the power series can be done directly:
Clear[seriesDSolve];
seriesDSolve[ode_, y_, iter : {x_, x0_, n_}, ics_: {}] :=
Module[{dorder, coeff},
dorder = Max[0, Cases[ode, Derivative[m_][y][x] :> m, Infinity]];
(coeff =
Fold[Flatten[{#1,
Solve[#2 == 0 /. #1,
Union@Cases[#2 /. #1, Derivative[m_][y][x0] /; m >= dorder,
Infinity]]}] &, ics,
Table[SeriesCoefficient[
ode /. Equal -> Subtract, {x, x0, k}], {k, 0, Max[n - dorder, 0]}]];
Series[y[x], iter] /. coeff) /; dorder > 0]
seriesDSolve[y'[x] == x^2 + (y[x])^2, y, {x, 0, 15}, {y[0] -> 1} ]
$$1 + x + x^2 + \frac{4\,x^3}{3} + \frac{7\,x^4}{6} + \frac{6\,x^5}{5} + \frac{37\,x^6}{30} + \frac{404\,x^7}{315} + \frac{369\,x^8}{315} + \frac{1961\,x^{10}}{1400} +$$
$$\frac{75092\,x^{11}}{51975} + \frac{1238759\,x^{12}}{831600} + \frac{9884\,x^{13}}{6435} + \frac{17121817\,x^{14}}{10810800} + \frac{115860954\,x^{15}}{70945875} + O[x]^{16}$$
Then we plot the true solution and its Adomian approximation with N=2 terms
sol = (y[x] /.
DSolve[{y'[x] == x^2 + (y[x])^2, y[0] == 1}, y[x], x])[[1]]

As we see from this graph, we can get a good approximation within a small interval [0,0.2]. To get a better approximation, however, we need to to add more terms. The solution has a singular point near x = 1, which we find by equating the denominator to zero in the explicit representation of the true solution:

FindRoot[(x (BesselJ[1/4, x^2/2] Gamma[1/4] - 2 BesselJ[-(1/4), x^2/2] Gamma[3/4])) == 0, {x, 1}]
{x -> 0.969811}
Therefore, the true solution blows up around the point 0.969811, but the power series approximation cannot detect the exact point no matter how many terms we use.

We show how the ADM works when we split the initial condition between several few first Adomian components. So we choose the initial approximation as the solution of the following initial value problem:
$u' = x^2 , \qquad u(0) = 1/3 .$
Its solution is
DSolve[{y'[x] == x^2 , y[0] == 1/3}, y[x], x]
$u_0 (x) = \frac{1}{3} \left( 1 + x^3 \right) .$
Since all Adomian polynomails were determined previously, we find the next two components in the Adomian series by solving the following initial value problems:
DSolve[{y'[x] == (1+x^3)^2 /9, y[0] == 1/3}, y[x], x] // Flatten
u[1][x_] = y[x] /. %
1/126 (42 + 14 x + 7 x^4 + 2 x^7)
u[0][x_]= (1+x^3)/3
A[1][x_] := Simplify[2*u[0][x]*u[1][x]]
DSolve[{y'[x] == A[1][x], y[0] == 1/3}, y[x], x] // Flatten
u[2][x_] = y[x] /. %
phi2[x_] = Simplify[u[0][x] + u[1][x] + u[2][x]]
1 + x/3 + x^2/27 + x^3/3 + x^4/9 + x^5/45 + x^7/63 + x^8/168 + ( 2 x^11)/2079
$\phi_2 (x) = u_0 + u_1 + u_2 = 1 + \frac{x}{3} + \frac{x^2}{27} + \frac{x^3}{3} + \frac{x^4}{9} + \frac{x^5}{45} + \frac{x^7}{63} + \frac{x^8}{168} + \frac{2\,x^{11}}{2079} .$
The next componet will be
A[2][x_] = Simplify[2*u[0][x]*u[2][x] + (u[1][x])^2]
DSolve[{y'[x] == A[2][x], y[0] == 0}, y[x], x] // Flatten
u[3][x_] = y[x] /. %
(1/3492720)(1164240 x + 388080 x^2 + 43120 x^3 + 194040 x^4 + 155232 x^5 + 30184 x^6 + 20790 x^8 + 9856 x^9 + 1855 x^12 + 208 x^15)
phi3[x_] = Simplify[phi2[x] + u[3][x]]
1 + (2 x)/3 + (4 x^2)/27 + (28 x^3)/81 + x^4/6 + x^5/15 + ( 7 x^6)/810 + x^7/63 + x^8/84 + (8 x^9)/2835 + (2 x^11)/2079 + ( 53 x^12)/99792 + (13 x^15)/218295
A[3][x_] = Simplify[2*u[0][x]*u[3][x] + 2*u[1][x]*u[2][x]]
DSolve[{y'[x] == A[3][x], y[0] == 0}, y[x], x] // Flatten
u[4][x_] = y[x] /. %
phi4[x_] = Simplify[phi3[x] + u[4][x]]
1 + (8 x)/9 + (10 x^2)/27 + (32 x^3)/81 + (83 x^4)/486 + ( 2 x^5)/15 + (14 x^6)/405 + (163 x^7)/8505 + x^8/56 + (8 x^9)/945 + ( 143 x^10)/113400 + (2 x^11)/2079 + (53 x^12)/49896 + ( 17 x^13)/57915 + (13 x^15)/218295 + (3613 x^16)/83825280 + ( 46 x^19)/12442815
\begin{align*} \phi_4 (x) &= u_0 + u_1 + u_2 + u_3 + u_4 \\ &= 1 + \frac{8\,x}{9} + \frac{10\,x^2}{27} + \frac{32\,x^3}{81} + \frac{83\,x^4}{486} + \frac{2\,x^5}{15} \\ & \quad + \frac{14\,x^6}{405} + \frac{163\,x^7}{8505} + \frac{x^8}{56} + \frac{8\,x^9}{945} + \cdots . \end{align*}
So we see that the standard ADM gives a better approximation than the preceding algorithm.

We perform the ADM approximation based on the accelerated Adomian polynomials. The initial Adomian polynomail is the same:
u[0][x_] = (1 + x/3)^3 DSolve[{y'[x] == (1 + x^3/3)^2, y[0] == 0}, y[x], x] // Flatten u[1][x_] = y[x] /. %
1/126 (126 x + 21 x^4 + 2 x^7)}
J[1][x_]=(u[0][x]+u[1][x])^2 - (u[0][x])^2 DSolve[{y'[x] == J[1][x], y[0] == 0}, y[x], x] // Flatten u[2][x_] = y[x] /. %
(1/5239080)(5239080 x^2 + 1746360 x^3 + 1047816 x^5 + 291060 x^6 + 93555 x^8 + 34650 x^9 + 5040 x^11 + 2310 x^12 + 88 x^15)
And the third component becomes
J[2][x_] = (u[0][x] + u[1][x] + u[2][x])^2 - (u[0][x] + u[1][x])^2 DSolve[{y'[x] == J[1][x], y[0] == 0}, y[x], x] // Flatten u[3][x_] = y[x] /. %
(1/5239080)(5239080 x^2 + 1746360 x^3 + 1047816 x^5 + 291060 x^6 + 93555 x^8 + 34650 x^9 + 5040 x^11 + 2310 x^12 + 88 x^15)
phi3[x_] = Simplify[u[0][x]+u[1][x]+u[2][x] + u[3][x]]
1 + x + 2 x^2 + x^3 + x^4/6 + (2 x^5)/5 + x^6/9 + x^7/63 + x^8/28 + ( 5 x^9)/378 + (4 x^11)/2079 + x^12/1134 + (2 x^15)/59535
J[3][x_] = (u[0][x] + u[1][x] + u[2][x] + u[3][x])^2 - (u[0][x] + u[1][x] + u[2][x])^2 DSolve[{y'[x] == J[1][x], y[0] == 0}, y[x], x] // Flatten u[4][x_] = y[x] /. %
phi4[x_] = Simplify[phi3[x] + u[4][x]]
1 + x + 3 x^2 + (4 x^3)/3 + x^4/6 + (3 x^5)/5 + x^6/6 + x^7/63 + ( 3 x^8)/56 + (5 x^9)/252 + (2 x^11)/693 + x^12/756 + x^15/19845
■

Example: Consider the initial value problem for the Riccati equation

$y' = x^2 - y^2 , \qquad y(0) =1.$
We can use Mathematica to find its solution.
DSolve[{y'[x] == x^2 - (y[x])^2, y[0] == 1}, y[x], x]
{{y -> Function[{x}, ((1/2 + I/ 2) ((1 + I) x^2 BesselJ[-(3/4), (I x^2)/2] Gamma[1/4] + I Sqrt[2] x^2 BesselJ[-(5/4), (I x^2)/2] Gamma[3/4] + Sqrt[2] BesselJ[-(1/4), (I x^2)/2] Gamma[3/4] - I Sqrt[2] x^2 BesselJ[3/4, (I x^2)/2] Gamma[3/4]))/(x (BesselJ[1/4, ( I x^2)/2] Gamma[1/4] + (1 + I) Sqrt[2] BesselJ[-(1/4), (I x^2)/2] Gamma[3/4]))]}}
This solution is not yet useful because it contains the imaginary unit. We need to transfer it to a real value expression.
$y(x) = x\,\frac{I_{-3/4} \left( \frac{x^2}{2} \right) \pi \left( \Gamma^2 \left( \frac{3}{4} \right) \sqrt{2} + \pi \right) - 2\,K_{3/4} \left( \frac{x^2}{2} \right) \Gamma^2 \left( \frac{3}{4} \right)}{\pi \left( \Gamma^2 \left( \frac{3}{4} \right) \sqrt{2} + \pi \right) I_{1/4} \left( \frac{x^2}{2} \right) + 2\,K_{1/4} \left( \frac{x^2}{2} \right) \Gamma^2 \left( \frac{3}{4} \right)}$
where Iν(z) is the modified Bessel function of the first kind and Kν(z) is the modified Bessel function of the second kind (alternatively, it is also called after William Thomson, 1st Baron Kelvin); here $$\Gamma (\nu ) = \int_0^{\infty} x^{\nu -1} e^{-x}\,{\text d}x$$ is the gamma function.

The denominator in the formula above is never zero, as its graph shows. Therefore, the solution to the given IVP exists for all positive x.

deno[x_] = Pi*(Sqrt[2]*(Gamma[3/4])^2 + Pi)* BesselI[1/4, x^2 /2] + 2*(Gamma[3/4])^2 * BesselK[1/4, x^2 /2]
Plot[deno[x], {x, 0, 1.5}, PlotStyle -> Thick, PlotRange -> {{0.01, 1.4}, {14, 25}}]

The Maclaurin power series expansion of the solution can be found with Mathematica:

AsymptoticDSolveValue[{y'[x] == x^2 - (y[x])^2, y[0] == 1}, y[x], {x, 0, 10}]
1 - x + x^2 - (2 x^3)/3 + (5 x^4)/6 - (4 x^5)/5 + (23 x^6)/30 - ( 236 x^7)/315 + (201 x^8)/280 - (218 x^9)/315 + (2803 x^10)/4200
$y(x) = 1-x + x^2 - \frac{2}{3}\,x^3 + \frac{5}{6}\,x^4 - \frac{4}{5}\, x^5 + \frac{23}{30}\,x^6 - \frac{236}{315}\, x^7 + \frac{201}{280}\, x^8 - \frac{218}{315}\,x^9 + \frac{2803}{4200}\,x^{10} - \frac{33412}{51975}\,x^{11} + \cdots .$
According to G. Adomian, we seek the solution in a way similar to the previous example:
$y(x) = \sum_{n\ge 0} u_n ,$
where the initial term u0 is exactly the same: $$u_0 = 1 + x^3 /3 .$$ Next, we find five Adomian polynomials using Mathematica
f[u_] := - u*u;
{x^2 - u[0]^2, -2 u[0] u[1], -u[1]^2 - 2 u[0] u[2], -2 u[1] u[2] - 2 u[0] u[3], -2 (u[2]^2/2 + u[1] u[3]) - 2 u[0] u[4], -2 (u[2] u[3] + u[1] u[4]) - 2 u[0] u[5]}
\begin{align*} A_0 &= - u_0^2 , \\ A_1 &= -2\,u_0 u_1 , \\ A_2 &= -u_1^2 -2\,u_0 u_2 , \\ A_3 &= -2\, u_1 u_2 - 2\,u_0 u_3 , \\ A_4 &= - u_2^2 - 2\,u_1 u_3 -2\,u_0 u_4 , \\ A_5 &= -2\,u_2 u_3 - 2\, u_1 u_4 -2\,u_0 u_5 . \end{align*}
We proceed in exactly the same way as in the previous example. Now we are ready to determine the first five terms in the Adomian series by solving the recursive initial value problems:
$\frac{{\text d} u_n}{{\text d}x} = A_{n-1} \left( u_0 , u_1 , \ldots , u_{n-1} \right) , \qquad u_n (0) =0 , \qquad n=1,2,\ldots .$
u[0][x_] = 1 + x^3 /3 ;
DSolve[{y'[x] == - (u[0][x])^2, y[0] == 0}, y[x], x] // Flatten ;
u[1][x_] = y[x] /. %
{1/126 (-126 x + 42 x^3 - 21 x^4 - 2 x^7)}
$u_1 (x) = \frac{1}{126} \left( -126 x + 42 x^3 - 21 x^4 - 2 x^7 \right) .$
The second term is
DSolve[{y'[x] == -2*u[0][x]*u[1][x], y[0] == 0}, y[x], x] // Flatten;
u[2][x_] = y[x] /. %
(83160 x^2 - 13860 x^4 + 16632 x^5 - 2640 x^7 + 1485 x^8 + 80 x^11)/83160
So
$u_2 (x) = \frac{1}{83160} \left( 83160 x^2 - 13860 x^4 + 16632 x^5 - 2640 x^7 + 1485 x^8 + 80 x^{11} \right) .$
Now we add the first two terms:
phi2[x_] = Simplify[u[0][x] + u[1][x] + u[2][x]]
1 - x + x^2 + (2 x^3)/3 - x^4/3 + x^5/5 - x^7/21 + x^8/56 + ( 2 x^11)/2079
Next approximation:
DSolve[{y'[x] == -2*u[0][x]*u[2][x] - u[1][x]*u[1][x], y[0] == 0}, y[x], x] // Flatten;
u[3][x_] = y[x] /. %
(1/3492720)(-3492720 x^3 + 698544 x^5 - 814968 x^6 - 55440 x^7 + 124740 x^8 - 88704 x^9 + 10080 x^11 - 5565 x^12 - 208 x^15)
phi3[x_] = Simplify[phi2[x] + u[3][x]]
1 - x + x^2 - (2 x^3)/3 - x^4/6 + x^5/5 - ( 7 x^6)/30 - x^7/63 + x^8/56 - (8 x^9)/315 + (2 x^11)/2079 - ( 53 x^12)/33264 - (13 x^15)/218295
Hence,
\begin{align*} \phi_3 (x) &= u_0 (x) + u_1 (x) + u_2 (x) +u_3 (x) = 1 - x + x^2 - \frac{2\,x^3}{3} \\ & \quad - \underbrace{\frac{x^4}{6} + \frac{x^5}{5} - \frac{7\, x^6}{30} - \frac{x^7}{63} + \frac{x^8}{56} - \frac{8\, x^9}{315} + \frac{2\, x^{11}}{2079} - \frac{53\, x^{12}}{33264} - \frac{13\, x^{15}}{218295}}_{\mbox{noise terms}} . \end{align*}
Then we plot the true solution along with Adomian's approximation:
exact[x_] = Simplify[y[x] /. DSolve[{y'[x] == x*x - y[x]*y[x], y[0] == 1}, y[x], x]];

As we see, the Adomian approximation with three terms ϕ3 matches the true solution up to power 3, and the rest are noise terms that contain x powers greater than 3. We eliminate some noise terms with the next iteration.
DSolve[{y'[x] == -2*u[1][x]*u[2][x] - 2*u[0][x]*u[3][x], y[0] == 0}, y[x], x] // Flatten;
u[4][x_] = y[x] /. %
(1/103524220800)(103524220800 x^4 + 27606458880 x^7 + 3524753232 x^10 + 273490560 x^13 + 13386165 x^16 + 382720 x^19)
phi4[x_] = Simplify[phi3[x] + u[4][x]]
1 - x + x^2 - (2 x^3)/3 + (5 x^4)/6 + x^5/5 - (7 x^6)/30 + ( 79 x^7)/315 + x^8/56 - (8 x^9)/315 + (143 x^10)/4200 + ( 2 x^11)/2079 - (53 x^12)/33264 + (17 x^13)/6435 - ( 13 x^15)/218295 + (3613 x^16)/27941760 + (46 x^19)/12442815
which gives the correct approximation up to the power 4 accompanied by noise terms of powers greater than 4. The next term in the Adomian series:
DSolve[{y'[x] == -2*u[1][x]*u[2][x] - 2*u[0][x]*u[3][x], y[0] == 0}, y[x], x] // Flatten;;
u[5][x_] = y[x] /. %
(1/18700822293753600)(-18700822293753600 x^5 - 5610246688126080 x^8 - 819274119535872 x^11 - 74948248743984 x^14 - 4598143925760 x^17 - 185499373059 x^20 - 4293552640 x^23)
phi5[x_] = Simplify[phi4[x] + u[5][x]]
1 - x + x^2 - (2 x^3)/3 + (5 x^4)/6 - (4 x^5)/5 - (7 x^6)/30 + ( 79 x^7)/315 - (79 x^8)/280 - (8 x^9)/315 + (143 x^10)/4200 - ( 2227 x^11)/51975 - (53 x^12)/33264 + (17 x^13)/6435 - ( 43327 x^14)/10810800 - (13 x^15)/218295 + (3613 x^16)/27941760 - ( 1318 x^17)/5360355 + (46 x^19)/12442815 - (7523 x^20)/758419200 - ( 15178 x^23)/66108676095
We compare our ADM approximations with the true solution and the approximation obtained with the aid of the classical Runge--Kutta of order 4 (abbreviated by RK4) with step size h = 0.1.
numerator[x_] = x*((Sqrt[2]*(Gamma[3/4])^2 + Pi)*BesselI[-3/4, x^2/2]*Pi - 2*(Gamma[3/4])^2*BesselK[3/4, x^2/2])
deno[x_] = Pi*(Sqrt[2]*(Gamma[3/4])^2 + Pi)*BesselI[1/4, x^2/2] + 2*(Gamma[3/4])^2*BesselK[1/4, x^2/2]
phi[x_] = numerator[x]/deno[x]
Clear[y]; h = 0.1;
f[x_, y_] := x^2 - y^2
y[0] = 1;
Do[k1 = h f[h n, y[n]];
k2 = h f[h (n + .5), y[n] + k1/2];
k3 = h f[h (n + .5), y[n] + k2/2];
k4 = h f[h (n + 1), y[n] + k3];
y[n + 1] = y[n] + 1/6*(k1 + 2 k2 + 2 k3 + k4), {n, 0, 4}]
y[1]
0.835786
x = 0.1 0.909409 0.909418 0.909408 0.90941
x = 0.2 0.835785 0.836052 0.835732 0.835732
x = 0.3 0.777238 0.779122 0.776672 0.777238

So we see that the 5-term ADM approximation is less accurate than the Runge--Kutta algorithm of order four; Also, the fomer requires more labor compared to the RK-4.    ■

Example: Consider the initial value problem for the autonomous equation

$y' = -\sqrt{1/4 + y^2} , \qquad y(0) =3.$
Its solution exists for all x and is given by
$y(x) = -\frac{1}{2}\,\sinh \left( x- \mbox{arcsinh}(6) \right) .$
Its series representation is
$y(x) = 3 - \frac{\sqrt{37}}{2}\,x + \frac{3\,x^2}{2} - \frac{\sqrt{37}}{12}\,x^3 + \frac{x^4}{8} - \frac{\sqrt{37}}{240}\,x^5 + \frac{x^6}{240} - \frac{\sqrt{37}}{10080}\,x^7 + \frac{x^8}{13440} - \frac{\sqrt{37}}{725760}\,x^9 + \frac{x^{10}}{1209600} - \frac{\sqrt{37}}{79833600}\,x^{11} + \cdots ,$
because
exact[x_] = -Sinh[x - ArcSinh[6]]/2 ;
Series[exact[x], {x, 0, 12}]
According to the ADM, we assume that the solution can be represented by a convergent series
$y(x) = \sum_{n\ge 0} u_n (x) ,$
where the initial approximation is taken to be u0 = 3 because the given differential equation is autonomous. Substituting the infinite series into the given differential equation, we obtain
$\texttt{D} \, \sum_{n\ge 0} u_n (x) = \sum_{n\ge 0} \texttt{D} \, u_n (x) = \sum_{n\ge 0} A_n , \qquad\mbox{with} \quad \texttt{D} = \frac{\text d}{{\text d}x} ,$
\begin{eqnarray*} A_0 &=& - \sqrt{\frac{1}{4} + u_0^2} = - \sqrt{\frac{1}{4} + 9} = -\frac{\sqrt{37}}{2} , \\ A_1 &=& - \frac{u_0 u_1}{\sqrt{1/4 + u_0^2}} = -\frac{6}{\sqrt{37}}\,u_1 , \\ A_2 &=& \frac{1}{2} \left[ \frac{u_0^2}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{1}{\sqrt{1/4 + u_0^2}} \right] u_1^2 - \frac{u_0 u_2}{\sqrt{1/4 + u_0^2}} = -\frac{6}{\sqrt{37}}\,u_2 - \frac{1}{37\sqrt{37}}\,u_1^2 , \\ A_3 &=& \frac{1}{2} \left[ \frac{u_0}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{u_0^3}{\left( 1/4 + u_0^2 \right)^{5/2}} \right] u_1^3 + \left[ \frac{u_0^2}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{1}{\sqrt{1/4 + u_0^2}} \right] u_1 u_2 - \frac{u_0 u_3}{\sqrt{1/4 + u_0^2}} \\ &=& -\frac{6}{\sqrt{37}}\,u_3 - \frac{2}{37\sqrt{37}}\,u_1 u_2 + \frac{12}{1369\sqrt{37}}\, u_1^3 , \\ A_4 &=& \frac{1}{8} \left[ \frac{5\,u_0^4}{\left( 1/4 + u_0^2 \right)^{7/2}} - \frac{6\, u_0^2}{\left( 1/4 + u_0^2 \right)^{5/2}} + \frac{1}{\left( 1/4 + u_0^2 \right)^{3/2}}\right] u_1^4 + \frac{3}{2} \left[ \frac{u_0}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{u_0^3}{\left( 1/4 + u_0^2 \right)^{5/2}} \right] u_1^2 u_2 + \frac{1}{2} \left[ \frac{u_0^2}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{1}{\left( 1/4 + u_0^2 \right)^{1/2}} \right] \left( u_2^2 + 2\, u_1 u_3 \right) - \frac{u_0 u_4}{\left( 1/4 + u_0^2 \right)^{1/2}} \\ &=& -\frac{6}{\sqrt{37}}\,u_4 - \frac{1}{37\sqrt{37}}\left[ u_2^2 + 2\,u_1 u_3 \right] + \frac{36}{1369\sqrt{37}}\, u_1^2 u_2 - \frac{143}{50653\sqrt{37}}\, u_1^4 , \\ A_5 &=& \frac{1}{120} \left[ - \frac{105\,u_0^5}{\left( 1/4 + u_0^2 \right)^{9/2}} + \frac{150\, u_0^3}{\left( 1/4 + u_0^2 \right)^{7/2}} - \frac{45\, u_0}{\left( 1/4 + u_0^2 \right)^{5/2}} \right] u_1^5 + \frac{1}{6} \left[ \frac{15\, u_0^4}{\left( 1/4 + u_0^2 \right)^{7/2}} - \frac{18\, u_0^2}{\left( 1/4 + u_0^2 \right)^{5/2}} + \frac{3}{\left( 1/4 + u_0^2 \right)^{3/2}}\right] u_1^3 u_2 + \frac{3}{2} \left[ \frac{u_0}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{u_0^3}{\left( 1/4 + u_0^2 \right)^{5/2}} \right] \left( u_1 u_2^2 + u_1^2 u_3 \right) + \left[ \frac{u_0^2}{\left( 1/4 + u_0^2 \right)^{3/2}} - \frac{1}{\left( 1/4 + u_0^2 \right)^{1/2}} \right] \left( u_2 u_3 + u_1 u_4 \right) - \frac{u_0u_5}{\left( 1/4 + u_0^2 \right)^{1/2}} \\ &=& - \frac{6}{\sqrt{37}}\,u_5 - \frac{2}{37\sqrt{37}} \left[ u_2 u_3 + u_1 u_4 \right] + \frac{36}{1369\sqrt{37}} \left[ u_1 u_2^2 + u_1^2 u_3 \right] - \frac{572}{50653\sqrt{37}}\, u_1^3 u_2 + \frac{1692}{1874161\sqrt{37}} \, u_1^5 , \end{eqnarray*}
and so on. Then we get a sequence of initial value problems:
$\frac{\text d}{{\text d}x} \, u_{n+1} = A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (0) =0, \quad n=0,1,2,\ldots .$
Consequently, we have
$\frac{\text d}{{\text d}x} \, u_{1} = -\frac{\sqrt{37}}{2} \qquad \Longrightarrow \qquad u_1 (x) = -\frac{\sqrt{37}}{2}\, x ,$
$\frac{\text d}{{\text d}x} \, u_{2} = \frac{6}{\sqrt{37}} \,\frac{\sqrt{37}}{2}\,x \qquad \Longrightarrow \qquad u_2 (x) = \frac{3}{2}\, x^2 ,$
Now we proceed similarly.
\begin{eqnarray*} \frac{\text d}{{\text d}x} \, u_{3} &=& A_2 \qquad \Longrightarrow \qquad u_3 (x) = - \frac{\sqrt{37}}{12}\, x^3 , \\ \frac{\text d}{{\text d}x} \, u_{4} &=& A_3 \qquad \Longrightarrow \qquad u_4 (x) = \frac{x^4}{8} , \\ \frac{\text d}{{\text d}x} \, u_{5} &=& A_4 \qquad \Longrightarrow \qquad u_5 (x) = - \frac{\sqrt{37}}{240} \, x^5 , \end{eqnarray*}
$\phi_5 (x) = u_0 + u_1 + u_2 + u_3 + u_4 + u_5 = 3 - \frac{\sqrt{37}}{2}\,x + \frac{3\,x^2}{2} - \frac{\sqrt{37}}{12}\,x^3 + \frac{x^4}{8} - \frac{\sqrt{37}}{240}\,x^5 .$
Finally we plot this approximation along with the true solution.
phi5[x_] = 3 - Sqrt[37]*x/2 + 3*x^2 /2 - Sqrt[37]*x^3 /12 + x^4 /8 - Sqrt[37]*x^5 /240 ;
Plot[{Callout[exact[x], "exact"], Callout[phi5[x], "ADM 5-approximation"]}, {x,0,3.0}, PlotStyle->Thick]
■

Since the ADM inverntion, many modifications and improvements were introduced in its applications, including several algorithms for evaluation of Adomian's polynomials and their generalizations. Instead of the standard ADM procedure, Abdul-Majid Wazwaz and some other researchers suggest splitting the inhomogeneous input between the first two (or more) terms in Adomian's series to minimize the influence of noise terms that are typical in nonhomogeneous equations. Another approach is based on splitting the initial condition between sequential initial value problems in order to increase the region of convergence for Adomian's decomposition solution. It represents the solution of the initial value problem
$y' = f(x,y) + g(x) , \qquad y(x_0 ) = y_0$
as the sum of components $$y(x) = u_0 + \sum_{n\ge 1} u_n (x) ,$$ where u0 and un, n ≥ 1, are solutions of the following initial value problems
$u'_0 = g(x) , \qquad u_0 (x_0 ) = c_0 ,$
and
$u'_n = A_{n-1} \left( u_0 , u_1 , \ldots , u_{n-1} \right) , \qquad u_n (x_0 ) = c_n , \quad n=1,2,\ldots ,$
where { cn }n≥0 is an arbitrary sequence of numbers adding to y0:
$y_0 = c_0 + c_1 + c_2 + \cdots = \sum_{n\ge 0} c_n .$
The following examples demonstrate implementation of the latter.

Example: Consider the initial value problem

$y' + y^2 =0, \qquad y(0) = 2.$
Its solution is
$y(x) = \frac{2}{2x+1} , \qquad x > 1/2.$
DSolve[{y'[x] == -y[x]*y[x], y[0] == 2}, y[x], x] // Flatten
{y[x] -> 2/(1 + 2 x)}
We decompose the solution and the nonlinearity into series:
$y(x) = \sum_{n\ge 0} u_n (x) \qquad\mbox{and} \qquad y^2 = \sum_{n\ge 0} A_n \left( u_0 , u_1 , \ldots , u_n \right) ,$
$A_n = \sum_{k=0}^n u_k u_{n-k} , \qquad n=0,1,2,\ldots .$
Upon introducing a parameter c < 2, we decompose the initial condition as
$y(0) = 2 = \left( 2-c \right) + c .$
Then the first two components u0 and u1 in the Adomian series are solutions of the following initial value problems
$u'_0 = 0 , \qquad u_0 (0) = 2-c \qquad \Longrightarrow \qquad u_0 (x) = 2-c;$
and
$u'_1 = -u_0^2 , \qquad u_1 (0) = c \qquad \Longrightarrow \qquad u_1 (x) =c - \left( 2-c \right)^2 x .$
DSolve[{y'[x] == -(2-c)^2, y[0] == c}, y[x], x] // Flatten
{y[x] -> c - 4 x + 4 c x - c^2 x}
Other components are solutions of the IVPs with homogeneous initial conditions
$u'_{n+1} = -A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (0) = 0 , \qquad n=1,2,\ldots .$
So \begin{align*} u_2 (x) &= \left( 2-c \right) \left[ -2cx + \left( 2-c \right)^2 x^2 \right] , \\ u_3 (x) &= -c^2 x + 3 \left( 2-c \right)^2 x^2 - \left( 2-c \right)^4 x^3 \\ u_4 (x) &= 6c^2 x^2 -3c^3 x^2 + \cdots . \end{align*}
u0=2-c
u1[x_] = c - 4 x + 4 c x - c^2 x
DSolve[{y'[x] == -(2*u0*u1[x]), y[0] == 0}, y[x], x] // Flatten
u[2][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[2][x] + u1[x]^2), y[0] == 0}, y[x], x] // Flatten
u[3][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[3][x] + 2*u1[x]*u[2][x]), y[0] == 0}, y[x], x] // Flatten
u[4][x_] = y[x] /. %
DSolve[{y'[x] == -(2*u0*u[4][x] + 2*u1[x]*u[3][x] + u[2][x]^2), y[0] == 0}, y[x], x] // Flatten
u[5][x_] = y[x] /. %
phi5[x_,c_]= Simplify[u0 + u1[x]+u[2][x]+u[3][x]+u[4][x]+u[5][x]]
The fifth degree polynomial approximation becomes
$\phi_5 (x,c) = u_0 (c) + u_1 (x,c) + u_2 (x,c) + u_3 (x,c) + u_4 (x,c) + u_5 (x,c) = \sum_{k=0}^5 u_k (x,c) .$
When c = 0, we get the classical Adomian series
$\phi (x) = y(x,0) = 2\sum_{n\ge 0} \left( -2x \right)^n ,$
which converges in the interval -1/2 < x < 1/2.

In order to see the effect of varying the value of c on the domain of convergence, we plot the curves of the explicit solution and fifth degree decomposition approximations for some values of c.

Plot[{2/(2*x + 1), phi5[x, 0], phi5[x, 1], phi5[x, 3/2]}, {x, -0.1, 1}, PlotStyle -> Thick, PlotLabels -> {"true", "c=0", "c=1", "c=3/2"}]
■

Example: Consider the initial value problem with negative power nonlinearity

$\frac{{\text d}y}{{\text d}x} = \frac{2}{y} , \qquad y(0) = 3 .$
Its explicit solution is
$y(x) = \sqrt{9+4x} = 3 + \frac{2}{3}\, x - \frac{2}{3^3}\,x^2 + \frac{4}{243}\,x^3 - \frac{10}{2187}\, x^4 + \frac{28\,x^5}{19 683} - \frac{28\,x^6}{59 049} +\cdots .$
DSolve[{y'[x] == 2/y[x], y[0] == 3}, y[x], x] // Flatten
Series[Sqrt[9 + 4*x], {x, 0, 10}
The Maclaurin series above for the true solution converges when |x| < 9/4 = 2.25. The coefficients of this power series are expressed through the binomial coefficients
$\sqrt{9+4x} = 3 \sum_{k\ge 0} \binom{1/2}{k} \left( \frac{4}{9}\,x \right)^k , \qquad \binom{1/2}{k} = \frac{\left( \frac{1}{2} \right)^{\underline{k}}}{k!} .$
By the ADM, we seek the solution to the given initial value problem in the form
$y(x) = u_0 (x) + \sum_{n\ge 1} u_n (x) ,$
where components un are determined recursively: \begin{align*} u'_0 (x) &= 0 , \qquad u_0 (0) = 3-c \qquad \Longrightarrow \qquad u_0 (x) = 3-c , \\ u'_1 (x) &= A_0 = \frac{2}{3-c} , \qquad u_1 (0) = c \qquad \Longrightarrow \qquad u_1 (x) = \\ u'_{n+1} (x) &= A_n \left( u_0 , u_1 , \ldots , u_n \right) , \qquad u_{n+1} (0) =0, \qquad n=1,2,3,\ldots . \end{align*}
u0[x_]= 3-c
DSolve[{y'[x] == 2/(3-c), y[0] == c}, y[x], x] // Flatten
u[1][x_] = (-3 c + c^2 - 2 x)/(-3 + c)
For determinantion of other components, we need the Adomian polynomials for the the nonlinearity 2/y: \begin{align*} A_1 &= -\frac{2\,u_1}{u_0^2} , \\ A_2 &= \frac{2\,u_1^2}{u_0^3} - \frac{2\,u_2}{u_0^2} , \\ A_3 &= \frac{4\,u_1 u_2}{u_0^3} - \frac{2\,u_1^3}{u_0^4} - \frac{2\,u_3}{u_0^2} , \\ A_4 &= \frac{2\,u_1^4}{u_0^5} - \frac{6\,u_1^2 u_2}{u_0^4} + \frac{2\,u_2^2 + 4\,u_1 u_3}{u_0^3} - \frac{2\,u_4}{u_0^2} , \\ A_5 &= \frac{8\,u_1^3 u_2}{u_0^5} - \frac{2\,u_1^5}{u_0^6} - \frac{6\,u_1 u_2^2 + 6\,u_1^2 u_3}{u_0^4} + \frac{4\,u_2 u_3 + 4\,u_1 u_4}{u_0^3} - \frac{2\,u_5}{u_0^2} . \end{align*}
f[y_] = 2/y
We use Mathematica to find other components
DSolve[{y'[x] == 2/(3 - c), y[0] == c}, y[x], x] // Flatten
u[1][x_] = y[x] /. %
(-3 c + c^2 - 2 x)/(-3 + c)
DSolve[{y'[x] == -2*u[1][x]/u0^2, y[0] == 0}, y[x], x] // Flatten
u[2][x_] = y[x] /. %
-((2 (-3 c x + c^2 x - x^2))/(-3 + c)^3)
DSolve[{y'[x] == 2*u[1][x]*u[1][x]/u0^3 - 2*u[2][x]/u0^2, y[0] == 0}, y[x], x] // Flattenn
u[3][x_] = y[x] /. %
-((2 (9 c^2 x - 6 c^3 x + c^4 x + 9 c x^2 - 3 c^2 x^2 + 2 x^3))/(-3 + c)^5)
DSolve[{y'[x] == 4*u[1][x]*u[2][x]/u0^3 - 2*u[1][x]^3 /u0^4 - 2*u[3][x]/u0^2, y[0] == 0}, y[x], x] // Flatten
u[4][x_] = y[x] /. %
-((2 (-27 c^3 x + 27 c^4 x - 9 c^5 x + c^6 x - 54 c^2 x^2 + 36 c^3 x^2 - 6 c^4 x^2 - 30 c x^3 + 10 c^2 x^3 - 5 x^4))/(-3 + c)^7)
DSolve[{y'[x] == 2*u[1][x]^4/u0^5 - 6*u[1][x]^2 *u[2][x] /u0^4 + 2*u[2][x]^2/u0^3 +4*u[1][x]*u[3][x]/u0^3 - 2*u[4][x]/u0^2, y[0] == 0}, y[x], x] // Flatten
u[5][x_] = y[x] /. %
-((2 (-27 c^3 x + 27 c^4 x - 9 c^5 x + c^6 x - 54 c^2 x^2 + 36 c^3 x^2 - 6 c^4 x^2 - 30 c x^3 + 10 c^2 x^3 - 5 x^4))/(-3 + c)^7)
We sum all these components to obtain fifth degree polynomial approximation
phi5[x_,c_]= u0 + u[1][x]+u[2][x]+ u[3][x]+u[4][x]+u[5][x]
Then we plot this function for different values of parameter c:
Plot[{Sqrt[9+4*x], phi5[x, 0], phi5[x, 1], phi5[x, 3/2]}, {x, -0, 1.5}, PlotStyle -> Thick, PlotLabels -> {"true", "c=0", "c=1", "c=3/2"}]
■
1. G. Adomian, K. Malakian, A comparison of the iterative method and Picard's successive approximations for deterministic and stochastic differential equations, Applied Mathematics and Computation, Vol. 8, No. 3, 187--204, 1981.
2. N. Bellomo, D. Sarafyan, On Adomian's decomposition method and some comparisons with Picard's iterative scheme, Journal of Mathematical Analysis and Applications, Vol. 123, No. 2, 389--400, 1987.
3. Duan, J.-S., Rach, R., Wang, Z., On the effective region of convergence of the decomposition series solution, Journal of Algorithms & Computational Technology, 2012, Vol. 7, No. 2, pp. 227--247.
4. Ray, S.S., New Approach for General Convergence of the Adomian Decomposition Method, World Applied Sciences Journal, 2014, Vol. 32, No. 11, pp. 2264--2268; doi: 10.5829/idosi.wasj.2014.32.11.1317