## Abstract

Exact solutions for nonlinear Arrhenius reaction–diffusion are constructed in *n*-dimensions. A single relationship between nonlinear diffusivity and the nonlinear reaction term leads to a non-classical Lie symmetry whose invariant solutions have a heat flux that is exponential in time (either growth or decay), and satisfying a linear Helmholtz equation in space. This construction also extends to heterogeneous diffusion wherein the nonlinear diffusivity factorizes to the product of a function of temperature and a function of position. Example solutions are given with applications to heat conduction in conjunction with either exothermic or endothermic reactions, and to soil–water flow in conjunction with water extraction by a web of plant roots.

## 1. Introduction

Assuming that volumetric heat capacity is constant and neglecting reagent consumption, the temperature of a reactive mixture satisfies a nonlinear reaction–diffusion equation
*Ω* is compact and connected, ∇ is the usual gradient operator and *t*_{2}>0. *R*(*θ*) is real-valued at any value of absolute temperature *D* must be positive [1]. One of the most important forms for the reaction function is the Arrhenius reaction term
*B* a positive constant. This follows from the Boltzmann–Gibbs equilibrium canonical distribution governing the probability of a particle overcoming the activation energy barrier *E* of a reaction. The parameter *B* may then be identified with *E*/*k*_{B}, where *k*_{B} is Boltzmann’s constant. The rate factor *ρCR*_{0}, *ρ* being density and *C* being specific heat at constant volume, is the amount of heat energy released per unit volume per unit time at very high temperature. This is usually described by a weak power-law dependence on temperature, *R*_{0}=*S*_{0}*θ*^{m}. For example, kinetic theory of a hard-sphere gas predicts *m*=0.5 [2]. Because the temperature-dependence is dominated by the Arrhenius factor, most extant models assume *m*=0. This will be assumed here, but a straightforward generalization of the following analysis can cover the case *m*≠0 with *ρC* also depending on *θ*.

For first-order reactions, *E* is the energy of bond dissociation into reactive molecular components such as free radicals. For example, for the exothermic decomposition of diethyl peroxide, measured reaction rates agree well with the Arrhenius law over a broad range of temperatures [3]. In particular, as a combustible mixture is controlled by extracting heat through the boundaries of the container, the Arrhenius reaction form is expected to remain appropriate.

A full Lie point symmetry classification of (1.1) was made by Dorodnitsyn *et al.* [4]. Various combinations of power-laws and exponential functions for *D*(*θ*) and *R*(*θ*) lead to an expansion of the Lie group of point symmetry transformations of (1.1) beyond the common Euclidean isometries in space, and translations in time. These, in turn, open possibilities of invariant solutions that may be obtained by reduction of variables. The invariant solutions indicate a wide range of possible dynamical behaviours, including stable similarity forms for the temperature with multi-peaks, extinguishing and blow-up in finite or infinite time that can occur even in solutions with compact support [4].

The non-analytic expression *θ*_{0} [6]. Wake & Bazley [7] extended the lower bound

Although the classical Lie point symmetry classification of (1.1), with *D*(*θ*) and *R*(*θ*) arbitrary, selects only constant or unbounded reaction functions, the complete *non-classical* symmetry classification of reaction–diffusion equations admits a broader range of possibilities [9–12], even admitting the Arrhenius reaction term [12]. Non-classical symmetries, in the sense of Bluman & Cole [13], leave invariant a system consisting of the original governing equation (1.1) plus the invariant surface condition that restricts the solution set to only invariant solutions. The concept of non-classical symmetry extends more generally to that of compatibility with an invariance condition, that may be of higher order than a point transformation or a contact transformation [14]. Non-classical symmetry analysis reveals that with a suitable nonlinear diffusivity function, after a change of variables, the Arrhenius reaction–diffusion equation admits separation of variables, resulting in a *linear* system. In the following §2, time-dependent radial solutions of this system are constructed in two- and three-dimensions. In the first example, the complicated nonlinear diffusivity function, as well as the temperature, are given exactly. With the same strongly increasing diffusivity, the exponentially heating solution is mathematically equivalent to that involving conduction through a finite domain with an endothermic reaction and exponential cooling. The strongly increasing diffusivity resembles soil–water diffusivity, with the sink term representing plant root absorption. In the case of ideal maximal cooling at the boundary, the nonlinear diffusivity is bounded, and it is the fixed point of a rapidly converging contraction map. The similarity form of the Kirchhoff variable is given exactly and it is asymptotically of the same form as the temperature. The diffusivity varies so little that the exact solution for the Kirchhoff variable closely approximates the temperature at all times. Particular attention is paid to the case of Newton cooling at the boundary, with non-zero Biot number.

In §3, it is shown that the same type of solution construction applies to a class of nonlinear reaction–diffusion equations with spatially varying diffusivity. The results of §2 are extended to allow for spatially variable diffusivity and a nonlinear sink term that applies to water transport in unsaturated soil with extraction by plant roots.

## 2. Non-classical reduction of nonlinear reaction–diffusion to the Helmholtz equation

A full non-classical symmetry classification of nonlinear diffusion–reaction equations in two spatial dimensions, was given in [12]. First, (1.1) is expressed in terms of the Kirchhoff variable [15],
*θ*_{0}≥0, then
*u*=0 corresponds to *θ*=*θ*_{0}.

In terms of the Kirchhoff variable, the reaction–diffusion equation is
*Q*(*u*)=*R*(*θ*) and *F*(*u*)=1/*D*(*θ*). The starting point of the reduction to the Helmholtz equation is the observation that equation (2.2) admits the non-classical group characterized by the reduction operator
*F* and *Q* are related by
*u*_{t}=*Au*. Then the symmetry reduction leads to

In terms of the original temperature variable *θ*, the relation (2.4) is
*R*(*θ*) from *D*(*θ*). Some basic combinations, (*D*(*θ*),*R*(*θ*)), with *R*(0)=0, are given in table 1. Note that the *R*(*θ*) function in case (d) of table 1 agrees asymptotically with the Arrhenius reaction term *R*_{0} *e*^{−B/θ} as inverse temperature approaches

However, in order to construct *D*(*θ*) from *R*(*θ*), it is necessary to solve a first-order nonlinear differential equation. From (2.4),
*u*(*θ*) the subject, then differentiating, this implies
*Φ* is simply the linear Helmholtz equation if *κ*=*K*^{2}>0, and the Laplace equation if *κ*=0. If *κ*=−*K*^{2}<0 with *R*>0, from (2.7), then *A* must be positive. Then the equation for *Φ* is the modified Helmholtz equation for which the radial solutions are unbounded either at the origin *r*=0 or at infinity.

### (a) An explicit solution with Arrhenius reaction: *κ*=0

In the case *κ*=0, equation (2.7) is linear homogeneous, allowing direct integration to obtain
*c*_{1} is an arbitrary constant that is positive (negative) if *R* is a strictly positive (negative) source (sink) term. However, the spatial dependence of *u*, given in (2.5), is then governed by Laplace’s equation. For example, the isotropic radial solutions *u*(**x**,*t*)=*u*(*r*,*t*) in two- or three-dimensions can only be a constant added to the singular point-source (sink) solutions with the source (sink) strength varying exponentially in time. The isotherms are specified exactly by the mapping that follows from equation (2.9):
*r*=*c*_{3}/(*Φ*−*c*_{2}), respectively.

For the case of the Arrhenius reaction (1.2), the solution (2.5) is valid when the diffusivity is exactly
*E*_{i} is the exponential integral. When *θ* is large, we can use the power series for the exponential integral
*γ* = Euler’s constant, to find
*θ* is small, we can use the asymptotic expression
*D*→0 as *θ*→0.

This nonlinear diffusivity is plotted in figure 1*b* alongside the Arrhenius reaction function, plotted in figure 1*a*.

### Example 2.1

A heat conduction problem with Arrhenius endothermic reaction term.

Even though the enthalpy of reaction is negative, an endothermic reaction will proceed naturally if the reaction products result in a lowering of the Gibbs free energy. An endothermic reaction may have a single activation energy, so that the rate of heat absorption is described by an Arrhenius function. Consider a spherical vessel containing the reagents of an endothermic reaction. *R*_{0} is negative. A small spherical decaying radioactive heat source at the centre supplies a quantity of heat Q. Then, because *u*>0 and *D*>0 when *θ*>0, (2.9)–(2.10) now imply *c*_{1}<0 and *A*<0. A solution
*u*.

### Example 2.2

Porous media flow with plant-root sink term.

The combination of exponentially growing nonlinear diffusivity, with bounded sink term, applies to water transport through a web of plant roots [16]. The domain of the problem is considered to be a cylindrical pile of soil of radius *r*=*r*_{1} with a vertical cylindrical injection well in the centre (*r*=*r*_{0}<*r*_{1}). The dependent variable *θ* now designates the water content above the plants’ wilting point, where the sink term approaches zero as roots fail to draw water. In this application, *R*_{0} is negative. Then, because *u*>0 and *D*>0 when *θ*>0, (2.9)–(2.10) now imply *c*_{1}<0 and *A*<0. A solution
*Q* per unit length of injection well into a large cylindrical soil mound that is exposed to the soil-controlled second stage of evaporation [17] at its outer boundary.

This solution applies analogously to a vertical current-carrying wire along the axis of a cylindrical region, supplying energy for an endothermic reaction.

### (b) Construction of bounded temperature solutions: *κ*>0

### Example 2.3

Heat conduction with exothermic reaction in a compact region.

With *κ*=*K*^{2}>0, *Φ*=*j*_{0}(*Kr*) in three-dimensions, *Φ*=*J*_{0}(*Kr*) in two-dimensions and *K*=λ_{1}/*r*_{1}, λ_{1} being the first zero of the spherical Bessel function *j*_{0} (λ_{1}=*π*), the standard Bessel function *J*_{0} (λ_{1}≈2.4048) or the cosine function (λ_{1}=*π*/2), respectively. Further, from the separation of variables in (2.5), that same solution satisfies other linear homogeneous boundary conditions such as
*Bi* constant. This resembles the linear condition for Newton cooling, in which *Bi* is the Biot number, after rescaling and non-dimensionalizing variables. The left-hand side is heat flux but the right-hand side is a multiple of the Kirchhoff variable rather than temperature. In the following, it is demonstrated that the nonlinear diffusivity has bounded variation and is close to being constant. This means that the above boundary condition is very close to that of Newtonian cooling to a very low-temperature environment. If the physical parameters *r*_{2},*D*(0) and *B*_{i} are prescribed, then the solution parameters are given by *A*=−*K*^{2}*D*(0), where *K* is the unique solution of −*Φ*′(*Kr*_{2})/*Φ*=*B*_{i}/*K*. *K* must lie between 0, where −*Φ*′(*Kr*_{2})/*Φ*=0 and λ_{1}/*r*_{2}, where −*Φ*′(*Kr*_{2})/*Φ* approaches infinity. With these homogeneous boundary conditions, there always exists an extinguishing solution in which the temperature approaches zero uniformly everywhere. A demonstration of the stability of this similarity solution is given in appendix A.

In the case *κ*>0, (2.7) is equivalent to the canonical form for an Abel equation of the second kind, via
*R*(*θ*) functions that produce solvable forms of equation (2.12). If *g*(*z*) a particular solution for an integrable form of (2.12), and *g*^{−1}(*w*) is the corresponding inverse function, an *R*(*θ*) function that produces (2.12) for *w*(*z*) is given implicitly according to

A good approximate analytic reconstruction of *D*(*θ*) may be obtained by applying a contraction map towards a fixed point. With *u*_{0}=0, the solution to (2.7) satisfying initial value *u*→0, *θ* → 0 must be a fixed point of the map
*A*| and *K* may be set to 1 by choosing dimensionless length and time coordinates *Kr* and |*A*|*t*.

In the case of the Arrhenius reaction term, we may set *R*(*θ*) to be *R*_{0} *e*^{−1/θ} by using a dimensionless temperature variable *θ*/*B*. Then the maximum value of *R*(*θ*)/*θ* is *R*_{0}/*e*.

From (2.14), using *D*(0)=1. From (2.14), it also follows that *D*(*θ*)>1 for *θ*>0. If *D*(*θ*) is bounded for *θ* sufficiently large, *D*(*θ*) must have a limit

Consider *f* on *f*−1 is non-negative
*D*_{n}=*M*^{n}*D*_{0}, *E*_{n}=*M*^{n}*E*_{0}, where *E*_{n+1}(*θ*) and

If *R*_{0}<*e*, then *R*_{0}=1 for which the above estimate of the contraction factor is 1/1.086. In practice, the convergence rate is better than this modest value suggests. In the case *D*_{0}(*θ*)=1. Using the standard notation *β*=1/*θ*, the iterative estimates are
*I* summation is bounded above by (*n*−2)!/*n*^{n−1}, which is of order *n*=5, *e*^{−5β} terms have a combined value of less than 0.01.

In appendix B, *D*(*θ*) is calculated accurately after obtaining exact series forms for *u*(*θ*). The simpler approximate diffusivity functions *D*_{j}(*θ*) are shown in figure 2 to agree well with the exact solution.

The diffusivity function has a single local maximum where its value is no higher than 45% above its mean value of 1. From (2.8), a stationary point of *D*(*θ*) may occur only where *D*(*θ*)=1+*R*′(*θ*)=1+*R*_{0}*θ*^{−2} *e*^{−1/θ}. This expression has a maximum value, implying

In effect, because at small-*t*, local disturbances extend by diffusion to a depth proportional to *θ*(*r*,*t*) asymptotically approaches the Kirchhoff variable *u*(*r*,*t*) which is identical to the temperature when *D*=1.

The construction of the nonlinear diffusivity function by a contraction map relied on the fact that with Arrhenius reaction, *βQ*(*β*) (=*R*(*θ*)/*θ*) has an upper bound less than 1. The same construction will apply to other reaction laws with this property, for example *m*>−1. This includes the case

### (c) The case *κ*<0

If *κ*=−*K*^{2}<0, from (2.7) with *R*(*θ*)>0, *A* must be positive. The diffusivity, *D*(*θ*), has a singularity where at leading order *D*′(*θ*)∼(*κ*^{2}/*A*)*D*^{3} and *D*∼(*θ*_{a}−*θ*)^{−1/2} at some freely chosen positive value *θ*=*θ*_{a}. In the modelling of unsaturated soil–water flow, the location of the singularity is chosen to be slightly greater than the volumetric water content at saturation [20].

Also the equation for *Φ* is the modified Helmholtz equation for which the radial solutions are unbounded either at the origin *r*=0 or at infinity. There is an exact solution for exothermic reaction–diffusion on the domain *r*>*r*_{0},
*K*_{0} is the modified Bessel function of the second kind. Although *u* is unbounded, the corresponding value of *θ* remains less than the singular value *θ*_{a}, owing to the singularity in *D*(*θ*) within the integrand in the Kirchhoff transformation *θ*→*u*. This singular nonlinear diffusivity does not occur in heat conduction but it may occur in models of unsaturated soil–water flow. However, in that application, positive distributed water sources are not common. It is more common to consider distributed sinks, owing to plant roots.

### Example 2.4

Porous media flow with plant-root sink term: *κ*<0.

As in example 2.2, we consider a distributed sink with *R*_{0}<0, owing to plant roots. Because of plant root extraction, *A*<0, so that a finite amount of water is supplied over all time through the inner surface at *r*=*r*_{0}. From (2.6), *D*>0,*A*<0 and *κ*<0 implies *R*<0. However in this case, from the solution (2.17), *θ* does not reach zero at any finite value of *r* but instead approaches zero at infinite distance. It is to be noted that the modified Helmholtz equation has previously been applied to *steady-state* solutions of unsaturated soil–water flow, by analogy to external problems of Helmholtz acoustic scattering [21]. Just as in the latter, we may construct analogous exterior solutions for boundaries of various shapes. However, at large values of *r*, it is well known that the solutions agree asymptotically with the isotropic ones, so the spherical or circular boundaries are canonical.

## 3. Incorporating heterogeneity

We now consider the case in which the medium (or substrate) is spatially heterogeneous. The appropriate governing reaction–diffusion equation can be written
*f* that is the amplification factor of the diffusivity. By rewriting this equation in terms of the Kirchhoff variable (2.1), we obtain
*Q*(*u*)=*R*(*θ*) and *F*(*u*)=1/*D*(*θ*), as before. Equation (3.1) admits the same simple non-classical symmetry (2.3), whenever *F*(*u*) and *Q*(*u*) are related by Equation (2.4). This means that *R*(*θ*) can be constructed from *D*(*θ*) using (2.6), or *D*(*θ*) could be constructed from *R*(*θ*) by solving (2.7). By again setting *u*(**x**,*t*)=*e*^{At}*Φ*(**x**), we derive the second-order linear equation for *Φ*(*x*),
*n*-dimensions. In the case where *A*<0, the solution for *θ* will be asymptotically similar to *u*. The relationship between *D*(*θ*) and *R*(*θ*) is the same as that described in the preceding section, and the solutions for *D*(*θ*) also hold in the heterogeneous examples described in the following.

### (a) The case *κ*=0

### Example 3.1

Porous media flow with a heterogeneous substrate: *κ*=0.

As an example, we now reconsider the porous medium plant-root extraction problem described in §2a, with the inclusion of spatial heterogeneity representing a changing scale of soil pores. In Miller self-similar soils, *f* is the geometric scale factor of the soil pores [22,23], with *f*(**r**_{0})=1 at some chosen reference point.

In this example, the problem is written in (one-dimensional) radially symmetric coordinates and the heterogeneity is described by *f*(**x**)=*f*(*r*). The differential equation for *Φ*(*r*) replacing *Φ*(**x**) in (3.2) is
*f*(*r*) to give
*c* constant. From (2.11a), we deduce that *A*<0. Having freely chosen *f*(*r*_{0}) to be 1, the solution, in terms of the Kirchhoff variable, satisfying boundary conditions (2.11), is

### (b) The case *κ*>0

### Example 3.2

Heat conduction in a heterogeneous medium.

Let *κ*=*K*^{2}. The differential equation for *Φ*(*r*) replacing *Φ*(**x**) in (3.2) is
*f*(*r*) any power-law of the cylindrical radius, exact solutions are readily available [18].

With heterogeneity described by the factor *f*(*r*)=*r*_{0}/*r*, (3.3) reduces to Airy’s equation, with general solution

When *f*(*r*) is of the form
*Φ*(*r*_{1})=0 by setting *Kr*_{0}>1, the general real solution is
*u*(*r*_{1})=0. However, *r*=0 is a singular point where thermal conductivity is either zero or infinite. The solution may be regarded as exterior to the surface of a hot wire at *r*=*r*_{0} from where a finite quantity of heat is supplied.

### (c) The case *κ*<0

### Example 3.3

Porous media flow with a heterogeneous substrate: *κ*<0.

With the same heterogeneity as in the previous model (3.4) but with *κ*=−*K*^{2}<0, the solution with *f*(*r*)=(*r*/*r*_{0})^{2} is simply
*f*(*r*)=*r*_{0}/*r*,

The parameters *c*_{i} and *K* may be chosen, so that the exterior solution satisfies the same boundary conditions on

## 4. Conclusion

Provided the nonlinear diffusivity and the nonlinear reaction term satisfy a single relationship, the Kirchhoff variable *u*, which is the integral of nonlinear diffusivity, admits solutions that are obtainable by separation of variables to a linear system, whose solution is an exponential in time, multiplying an arbitrary solution of the Helmholtz, modified Helmholtz or Laplace equation in space. The heat flux is merely −∇*u*, which is given explicitly everywhere in these solutions.

If the nonlinear diffusivity function is specified, then the compatible reaction function can be constructed directly by integration. If the nonlinear reaction function is specified, then the diffusivity function is the solution of a differential equation that is equivalent to an Abel equation if *u* satisfies a Helmholtz equation, or to a separable equation if *u* satisfies Laplace’s equation. To the best of our knowledge, this provides the only known exact closed-form solutions in two- and three-dimensions, of nonlinear reaction–diffusion equations with the classical Arrhenius reaction term. Even when *D*(*θ*) satisfies an Abel equation that is not of known solvable type, in some circumstances, it may be specified to arbitrary precision using exact series expansions.

This construction is valid in any natural number *n* of dimensions, and it generalizes also to heterogeneous extensions of the Helmholtz factor equation. Applications are given for radial solutions of exothermic reactions with nonlinear heat conduction, endothermic reactions with nonlinear heat conduction, and water flow from a supply well into a cylindrical soil mound with soil-limited evaporation at the outer boundary. The logistic shape of the negative Arrhenius reaction term resembles the behaviour of distributed plant roots that have a maximum and minimum value of water uptake rate near saturation and wilting point, respectively.

As is well known from acoustic scattering theory, exterior solutions of the Helmholtz equation typically asymptotically approach radial symmetric solutions at large distances from the scattering surface. Hence, the radial solutions illustrated here are in a sense, canonical. The solution method, used here, involves a free function of the Helmholtz equation, so it would not be difficult to use known non-radial solutions. The special solutions presented here could at least be used as bench tests for more broadly applicable approximate solution methods.

This approach will lead to ongoing investigations of other nonlinear partial differential equations with more than one free function, that may admit special non-classical symmetry reductions.

## Competing interests

We declare we have no competing interests.

## Funding

We received no funding for this study.

## Disclaimer

P.B. conceived the general approach, provided the solution in the case *K*=0, provided recursive approximations for diffusivity and directed the project. B.H.B.-H. provided the solutions for heterogeneous models. D.T. provided series constructions for nonlinear diffusivity. All authors contributed to the writing, and all agree on the current version.

## Acknowledgements

We thank Dr Joanna Goard, Prof. Joel Moitsheki and Dr Edoardo Daly for useful discussions on the symmetry reductions and plant-root absorption.

## Appendix A. Stability of similarity solution

It is shown here that the extinguishing radial similarity solution of the Arrhenius reaction–diffusion equation given in §2b, is exponentially stable to small perturbations. Change variables to a set of canonical variables of the non-classical symmetry, namely *r*,*t* and *v*=*u* *e*^{−At}. In terms of these variables, the reaction–diffusion equation may be expressed
*v*(*r*_{1})=0. The similarity solution is a pseudo-steady state, *v*_{s}=*Φ*(*r*)=*A*_{0}*J*_{0}(λ_{1}*r*/*r*_{1}), corresponding to exponential decrease of the Kirchhoff variable *u*_{s}=*Φ*(*r*) *e*^{−|A|t} and satisfying *Φ*′(0)=0 and *Φ*(*r*_{1})=0. Now consider a perturbed solution, in plane polar coordinates
*ϵ* perturbation *w* satisfies the same homogeneous boundary conditions, plus 2*π*-periodicity with respect to *ϕ*, plus
*θ*, this is a linear equation with no squared derivative terms, allowing recourse to comparison theory of reaction–diffusion equations [24]. Note that 1/*F*(*u*)=*D*(*θ*) and that from §2b, *D*_{m}≥*D*≥*D*(0)=−*A*/*K*^{2}. Hence by comparison, |*w*(*r*,*t*)| must decay at least as fast as the solution to the linear initial-boundary problem with smaller diffusivity and larger reaction term,
*q*_{τ}=∇^{2}*q* with *q* satisfying the same initial and boundary conditions as *Q*. The solution *q*(*r*,*ϕ*,*t*) may be expanded as a standard Fourier–Bessel series [25]. Without loss of generality, we neglect the first component in the series expansion of *w*_{0}, which is a multiple of *J*_{0}(λ_{1}*r*/*r*_{1}), and which takes the form of the assumed unperturbed similarity solution for *v*. Each of the other terms in the series for *Q* is of the form
_{n,m} is the *m*th zero of Bessel function *J*_{n}, whereas *A*_{n,m} and *B*_{n,m} are arbitrary real coefficients. Note that beyond λ_{0,1}, the next smallest root is λ_{1,1}. From the estimate (2.16), each of the terms (A 1) is decreasing exponentially in time provided

## Appendix B. Evaluating *D*(*θ*) using asymptotic series

The reaction term *κ*=*K*^{2}, the parameters *K*, |*A*| and *B* in equation (2.7) may be set to 1 by adopting appropriate dimensionless variables. By successive isolation of leading-order terms, the asymptotic behaviour of *u*(*θ*) near *θ*=0 can be shown to be
*u*_{0}=0 but extension to *u*_{0}>0 should be straightforward. Adopting *q*_{r,−1}=0, substitution into (2.7) shows that
*q*_{r,m}} determined iteratively as above, the series (A 4) is divergent, but can still be used to give accurate specifications of *u*(*θ*) for *θ* sufficiently small. Partial sums of (A 4) may be produced by evaluating both *r* and *m* sums up to some maximum integer. Performance of the partial sums is improved by substituting the known asymptotic form

A Taylor series expansion for *u*(*θ*) centred on a non-singular point _{0}=*u*(*θ*_{0}) is assumed known, and the above sum evaluates to zero for *i*=0. The asymptotic series form of *u*(*θ*) as

We can use the above series to evaluate *u* for *R*_{0}=1 and *A*<0. Using (A 4) with *θ* < 0.9, and matches the accuracy of the known value at *θ*_{0}=10, with *u*(10)=11.79313028084656, extends the domain of the solution to *θ*=19.5, and is compatible with the known value *D*(*θ*)=*du*/*dθ*. To obtain a solution of greater accuracy, we can begin with evaluation at a point *θ*_{*}≤19.5 as done above.

- Received August 20, 2015.
- Accepted October 27, 2015.

- © 2015 The Author(s)