# Asymptotically-consistent analytical solutions for the non-Newtonian
Sakiadis boundary layer

Nastaran Naghshineh<EMAIL_ADDRESS>School of Mathematical Sciences, Rochester
Institute of Technology, Rochester, NY, 14623, USA  Department of Sciences and
Liberal Arts, Rochester Institute of Technology-Dubai, Dubai, 341055, UAE
Nathaniel S. Barlow  School of Mathematical Sciences, Rochester Institute of
Technology, Rochester, NY, 14623, USA  Mohamed A. Samaha  School of
Mathematical Sciences, Rochester Institute of Technology, Rochester, NY,
14623, USA  Department of Mechanical and Industrial Engineering, Rochester
Institute of Technology-Dubai, Dubai, 341055, UAE  Steven J. Weinstein  School
of Mathematical Sciences, Rochester Institute of Technology, Rochester, NY,
14623, USA  Department of Chemical Engineering, Rochester Institute of
Technology, Rochester, NY, 14623, USA

###### Abstract

The Sakiadis boundary layer induced by a moving wall in a semi-infinite fluid
domain is a fundamental laminar flow field relevant to high speed coating
processes. This work provides an analytical solution to the boundary layer
problem for Ostwald-de Waele power law fluids via a power series expansion,
and extends the approach taken for Newtonian fluids ["On the use of
asymptotically motivated gauge functions to obtain convergent series solutions
to nonlinear ODEs", IMA J. of Appl. Math., (2023)] in which variable
substitutions (which naturally determine the gauge function in the power
series) are chosen to be consistent with the large distance behavior away from
the wall. Contrary to prior literature, the asymptotic behavior dictates that
a solution only exists in the range of power law exponents, $\alpha$, lying in
the range $0.5<\alpha\leq 1$. An analytical solution is obtained in the range
of approximately $0.74\leq\alpha<1$, using a convergent power series with an
asymptotically motivated gauge function. For power laws corresponding to
$0.5<\alpha<0.74$, the gauge function becomes ill-defined over the full
domain, and an approximate analytical solution is obtained using the method of
asymptotic approximants ["On the summation of divergent, truncated, and
underspecified power series via asymptotic approximants", Q. J. Mech. Appl.
Math., (2017)]. The approximant requires knowledge of two physical constants,
which we compute a priori using a numerical shooting method on a finite
domain. The utility of the power series solution is that it can be solved on
the entire semi-infinite domain and–in contrast to a numerical solution–does
not require a finite domain length approximation and subsequent domain length
refinement.

###### pacs:

##  I Introduction

The Sakiadis boundary layer[1] is a fundamental flow field in processes where
laminar liquid films are coated onto moving substrates [2]. One of its key
physical implications is in the area of high speed curtain coating, where the
boundary layer length is essential to the mechanism of hydrodynamic assist
that can suppress air entrainment [3]. In particular, its length determines
where the wetting line is located with respect to the main body of the curtain
flow. Depending on the relative speed of the substrate and curtain flow at its
bottom, the wetting line can lie directly underneath the curtain or can be
dragged forward (lower curtain flow and higher substrate speeds) or retarded
backward (high curtain flow and lower substrate speed). This wetting line
location determines whether the stagnation pressure from a tall liquid curtain
is sufficient to suppress the creation of an unstable air-bearing that leads
to uneven and bubble-laden coatings. The highest coating speeds occur when the
wetting line is located directly underneath the centerline of the curtain
itself. The taller the liquid curtain, the faster the ultimate coating speed,
provided that the wetting line location–again dictated by the Sakiadis
boundary layer–is optimally controlled [3]. In addition to its relevance to
coating, the Sakiadis boundary layer forms the basis for many studies
including moving elastic sheets involving various modes of heat and mass
transport[4, 5, 6, 7, 8, 9, 10]; in the past 5 years, the original Sakiadis
paper [1] has been cited over 400 times which demonstrates its continued
fundamental importance.

Figure 1 shows the configuration of the Sakiadis boundary layer problem with
the $x$-$y$ coordinate system as indicated; the fluid flow is assumed to be
invariant with the direction oriented out of the figure. Here, a flat wall is
moving with velocity, $u_{w}$, through an otherwise stationary generalized
Newtonian incompressible fluid of density, $\rho$, and viscosity, $\mu$.

Figure 1: Schematic of the Sakiadis boundary layer flow. Slowly moving or
otherwise stationary fluid with characteristic velocity scale $S$ is in
contact with a fast moving wall having speed of $u_{w}$ as indicated, where
$S<<u_{w}$. As a result, the velocity $u\to 0$ as $y\to\infty$ in the boundary
layer approximation to the flow equations. The boundary layer envelope is
defined in this paper as the locus of points for which $u/u_{w}=0.1$. In the
original papers of Sakiadis [1] and Fox et al. [11], fluid enters the domain
at point $\mathbb{A}$ through a slit, where the streamline is redrawn to be
vertical and coincident with a wall. In coating applications, the streamline
often aligns with an interface where point $\mathbb{A}$ is a moving contact
line [2].

The governing equations embody conservation of mass and momentum through the
two dimensional incompressible steady continuity and Navier-Stokes equations.
At high substrate speeds, velocity in the $x$-direction, $u$, is much larger
than that in the $y$-direction, $v$, and velocity gradients in the
$y$-direction dominate viscous forces in the boundary layer. These
assumptions, which result in a small slope of fluid points, lead directly to
Prandtl’s boundary layer equations—a nonlinear partial differential equation
(PDE) system—that apply to the Sakiadis flow. For fluids having a viscosity
that satisfies the Ostwald-de Waele power law dependence, the Sakiadis
boundary layer equations are expressed as [11]

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0,$ (1a)
$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial
y}=\frac{1}{\rho}\frac{\partial\tau_{xy}}{\partial y},$ (1b)
$\tau_{xy}=\mu\left(\frac{\partial u}{\partial
y}\right)~{},~{}\mu=K\left(-\frac{\partial u}{\partial y}\right)^{\alpha-1},$
(1c) $u=u_{w},~{}v=0~{}\textrm{at}\ y=0~{};~{}u\to 0~{}\textrm{as}\
y\to\infty.$ (1d)

In (1dc), $\tau_{xy}$ is the shear stress in the fluid, $K>0$ is the
consistency coefficient, $\alpha$ is the power law exponent, and $\mu$ is the
strain-rate dependent viscosity. Note that the rate of strain in the boundary
layer approximation, $\partial u/\partial y$, is negative in the flow domain;
thus, the magnitude of the rate of strain invokes a negative sign as indicated
in the viscosity expression in (1dc). Through the use of the stream function,
$\psi$, that satisfies the continuity equation (1a) (i.e.,
$u=\partial\psi/\partial y$ and $v=-\partial\psi/\partial x$), Fox et al. [11]
define similarity variables given as

$\eta=y\left(\frac{\rho u_{w}^{2-\alpha}}{Kx}\right)^{1/(\alpha+1)},$ (2a)
$\psi=\left(\frac{y~{}u_{w}}{\eta}\right)f(\eta),$ (2b)

and thus $u$ and $v$ are expressed as

$u=u_{w}\frac{df}{d\eta}~{},~{}~{}v=\frac{1}{1+\alpha}\left(\frac{Ku_{w}^{2\alpha-1}}{\rho
x\alpha}\right)^{\frac{1}{\alpha+1}}\left(\eta\frac{df}{d\eta}-f\right).$ (3)

Upon substitution of (2b) into the system (1d) and after rearrangement, Fox et
al. [11] obtain the non-Newtonian Sakiadis boundary layer problem given as

$\alpha(\alpha+1)\frac{d^{3}f}{d\eta^{3}}-(-\frac{d^{2}f}{d\eta^{2}})^{(2-\alpha)}f=0,~{}~{}~{}0\leq\eta<\infty,$
(4a) $f=0~{}\textrm{at}\ \eta=0,$ (4b) $\frac{df}{d\eta}=1~{}\textrm{at}\
\eta=0,$ (4c) $\frac{df}{d\eta}=0~{}\textrm{as}\ \eta\to\infty.$ (4d)

Note that for a Newtonian fluid ($\alpha=1$), the equation system (4d) is
identical to that of Blasius [12] except that the location of conditions (4c)
and (4d) are reversed[1]. Nevertheless, the solution of (4d) is not simply
related to that of Blasuis by translation, owing to its nonlinear governing
equation (4a).

Analytical solutions to the Newtonian Sakiadis problem have been examined by
Barlow et al. [13], who show that a power series solution about $\eta=0$ has a
finite radius of convergence and cannot bridge the physical domain
$\eta\in[0,\infty)$. Naghshineh et al. [14] provide a convergent power series
solution for the Newtonian problem in terms of exponential gauge functions
consistent with the asymptotic behavior of the solution away from the wall;
this behavior is given as

$f(\eta)\sim
C+\tilde{a}_{1}e^{-C\eta/2}+\tilde{a}_{2}\left(e^{-C\eta/2}\right)^{2}+O\left(\left(e^{-C\eta/2}\right)^{3}\right)~{}\textrm{as}\
\eta\to\infty,$ (5)

where $C>0$ (this is consistent with the wall motion inducing a net volumetric
flow in the positive $x$-direction in Fig. 1), $\tilde{a}_{1}$ and
$\tilde{a}_{2}$ are asymptotic constants, and by inspection,

$\lim_{\eta\rightarrow\infty}f(\eta)\equiv C.$ (6)

A Taylor series solution is given as

$g(\omega)=\sum_{n=0}^{\infty}\tilde{a}_{n}\omega^{n},$ (7a)
$\omega(\eta)=e^{\displaystyle-C\eta/2},$ (7b) $f(\eta)\equiv
g(\omega(\eta)),$ (7c)

where all coefficients $\tilde{a}_{n}$ are provided by Naghshineh et al. [14]
As written, the solution (7c) uses $\eta=\infty$ as an expansion point; a
slightly faster converging Taylor series expansion (in this transformed gauge
function111A gauge function is the usual independent variable in a series
expansion. For example, in $e^{x}=\sum(x^{n}/n!)$, $x$ is the guage function.
If we write $\textrm{exp}(e^{x})=\sum(e^{nx}/n!)$, $e^{x}$ is the gauge
function [34, 35]) may be developed about $\omega=1$ (i.e. $\eta=0$)[14]. Note
that all constants, including the constant $C$, may be determined via an
algorithm independent of any numerical information. The reader is referred to
Naghshineh et al. [14] for a review of literature relevant to the Newtonian
Sakiadis problem.

The objective of this work is to extend the approach of Naghshineh et al. to
obtain an analytical solution of the Sakiadis boundary layer for Ostwald-de
Waele power law fluids, i.e. the solution of the system (4d) for $\alpha\neq
1$. This is especially relevant, as many fluids used in thin-film coating
exhibit shear thinning character [2], for which $\alpha\in(0,1)$. In this
parameter range, note that the power law model (1dc) is deficient in that it
limits to an infinite viscosity, $\mu$, in (1dc), as the rate of strain
approaches zero ($\partial u/\partial y\to 0$); nevertheless, as the shear
stress, $\tau_{xy}$, remains finite in this limit, reasonable predictions may
still be made, such as in pipe or slot flows[16]. The power-law model is often
used to describe flows for shear thinning behaviour due to its simplicity.
However, in this paper, we demonstrate that there is a restricted range of
$\alpha$ values ($0.5<\alpha\leq 1$) for which a solution to the power-law
non-Newtonian Sakiadis problem given by (4d) exists. This is a mathematical
restriction of the power-law itself, which is an approximation to the true
behavior of a shear thinning fluids [17, 18]. For fluids that exhibit a power-
law dependence satisfying $\alpha<0.5$ or $\alpha>1$ over a region of shear
rate, calculations must be done with more sophisticated viscosity dependences
such as the Carreau Model [17, 18]. When these models are used in place of
$\mu$ in (1dc), a similarity variable cannot be identified, and the full PDE
system governing boundary layer flows must be solved [19].

This paper is organized as follows. In Sec. II, we first examine the solution
to the system (4d) via a standard power series expansion about $\eta=0$, and
find that it is divergent, as was found for Newtonian fluids. In Sec. III, we
then consider the asymptotic solution of the system (4d) as $\eta\to\infty$,
and use it to motivate a Taylor series expansion in terms of an alternative
gauge function, as was done for the Newtonian solution (7c) discussed above.
By judiciously choosing the location of the expansion point, we are able to
obtain a convergent series solution for power law exponents lying in the range
of approximately $0.74\leq\alpha<1$. For power laws corresponding to
$\alpha<0.74$, the gauge function becomes ill-defined over the full domain. As
a result, in Sec. IV an accurate approximate solution is obtained using the
method of asymptotic approximants [13]. Here, the two necessary constants are
determined a priori using a numerical shooting method. Note that this approach
is taken in the solution of the Falkner-Skan boundary layer equations in prior
works [20, 21]. The utility of the analytical forms are demonstrated in Sec.
V, by the ease with which streamlines and the velocity field may be extracted.
Concluding comments are provided in Sec. VI. Formulae used to manipulate the
nonlinear series expansion in this study are provided in Appendix A. The
shooting algorithm used to solve the non-Newtonian Sakiadis problem (4d)
numerically is provided in Appendix B, and relevant constants for the
presented power series and approximant are provided in Appendix C. Appendix D
includes the algorithm used to predict the same constants via the convergent
power series solution itself.

##  II Power series solution

For this section and the next, we solve the ODE (4d), which arises after
similarity transform. As such, the physics of the original system (1d) is
obscured in the mathematical solution. In Sec. V, we demonstrate the ease with
which streamlines and velocity fields may be extracted from the solution we
provide, and in doing so provide solution results in the physical domain.

A power series solution to the ODE (4d) can be obtained through standard means
using JCP Miller’s formula [22] and Cauchy’s product rule [23] (see Appendixes
A.1 and A.2, respectively) to re-order nonlinear terms in powers of $\eta$;
the series expansion is

$f=\sum_{n=0}^{\infty}a_{n}\eta^{n},~{}~{}|\eta|<\eta_{s}(\alpha),$ (8a)
$a_{n+3}=\frac{\displaystyle\sum_{j=0}^{n}b_{j}~{}a_{n-j}}{\alpha(\alpha+1)(n+3)(n+2)(n+1)},~{}~{}n\geq
0,$ (8b) $b_{n>0}=\frac{1}{2na_{2}}\sum_{j=1}^{n}(3j-\alpha
j-n)(j+2)(j+1)a_{j+2}b_{n-j}~{}~{},~{}~{}b_{0}=(-2a_{2})^{2-\alpha},$ (8c)
$a_{0}=0~{}~{},~{}~{}a_{1}=1~{}~{},~{}~{}\textrm{and}\ a_{2}=\kappa/2,$ (8d)

where $\eta_{s}(\alpha)$ is a finite radius of convergence. In (8dd), the
quantity $\kappa$ is directly related to the wall shear stress in the boundary
layer flow, typically referred to as the "wall shear" parameter [24, 25], and
is defined as

$\kappa=f^{\prime\prime}(0).$ (9)

The quantity $\kappa$ in (9) is a function of $\alpha$, and is not known a
priori; it is typically determined numerically. Alternatively, $\kappa$ can be
calculated algorithmically as shown in Sec. III.3 of this paper as an
extension of the technique developed for the Newtonian Sakiadis problem[14].

Figure 2 provides a comparison between the power series solution (8d) and the
numerical solution to (4d) for $\alpha=0.8$. The numerical solution is
obtained using a shooting method (see Appendix B) to recast system (4d) as a
boundary value problem on a finite domain length $L$, where the condition (4d)
is replaced with $df/d\eta=0$ at $\eta=L$; the length $L$ is chosen such that
doubling its size leads to difference in the predictions of $\kappa$ and
$C=f(L)$ of $O(10^{-15})$ and $O(10^{-10})$ when $\alpha=0.8$. The constants
$\kappa$ and $C$ are defined in (9) and (6), respectively. As shown in the
Fig. 2, the power series solution (8d) diverges within the physical domain.
The rightmost vertical line marked by arrow $A$ shows the radius of
convergence of the power series solution (8d) given by $\eta_{s}(0.8)\approx
3.09$ , which is confirmed via a numerical root test 222Some of the
coefficients $a_{n}$ in (8d) are zero; hence, we use root test instead of
ratio test., as shown in Fig. 3.

Figure 2: The solution to (4d) is shown for $\alpha=0.8$. The numerical
solution (Appendix B) with $L=11000$ (black dots) is compared against the
$N$-term truncations of the power series solution (8d) (dashed curves),
transformed series solution (expansion about $\eta=\infty$) (16h) (dotted
curves), and the transformed series solution (expansion about $\eta=0$) (18j)
(solid curves) for $N=10$ and $N=25$. In regions of the plot where a given
dashed or dotted curve is not clearly seen, the curves agree with the
numerical results. The rightmost vertical line (marked by arrow A) shows the
radius of convergence $\eta_{s}\approx 3.09$ of the standard series solution
(8d), and the leftmost vertical line (marked by arrow B) shows the radius of
convergence $\eta_{s}\approx 0.95$ of the transformed series solution (16h).
For $\alpha=0.8$, the numerically obtained values of the constants $\kappa$,
$C$, and $E$ (defined in Sec.III) used in producing the figure are given in
Appendix C.2. Figure 3: Root test for (8d) is shown for $\alpha=0.8$,
indicating a radius of convergence ($y$-axis) of $\eta_{s}\approx 3.09$. This
is consistent with the divergent behavior observed in Fig. 2.

##  III Asymptotically motivated Gauge Function and Expansions

###  III.1 Asymptotic Behavior as $\eta\rightarrow\infty$

Similar to the approach taken for Newtonian fluids [14], we use the method of
dominant balance [27] to determine the asymptotic behavior of $f$ as
$\eta\to\infty$. This behavior motivates the use of a gauge function that
ultimately leads to a convergent series expansion. To proceed, we write the
solution of (4da) as

$f\sim C+h(\eta),~{}\textrm{with}\ h\to 0~{}\textrm{as}\ \eta\to\infty,$ (10a)
where $C$ is the asymptotic constant described in (6), and $h(\eta)$ is a
function to be determined. The form (10a) is substituted in (4a) to obtain
$\alpha(\alpha+1)h^{\prime\prime\prime}\sim(C+h)(-h^{\prime\prime})^{(2-\alpha)}~{}\textrm{as}\
\eta\to\infty,$ (10b) where the primes denote derivatives of $h$ with respect
to $\eta$. Equation (10b) may be simplified by noting that $h$ is subdominant
to $C$ as $\eta\to\infty$, and thus
$\alpha(\alpha+1)h^{\prime\prime\prime}\sim
C(-h^{\prime\prime})^{(2-\alpha)}~{}\textrm{as}\ \eta\to\infty.$ (10c) The
above equation can be integrated once to obtain
$h^{\prime\prime}\sim-\left[E+\frac{C(1-\alpha)}{\alpha(\alpha+1)}\eta\right]^{\frac{1}{\alpha-1}}~{}\textrm{as}\
\eta\to\infty,$ (10d) where $E$ is the constant of integration to be
determined. Integrating (10d) twice, and applying the boundary condition (4d),
the solution of (10d) is
$h\sim\displaystyle\frac{\alpha(\alpha+1)^{2}E^{\frac{2\alpha-1}{\alpha-1}}}{-C^{2}(2\alpha-1)}\left[1+\frac{C(1-\alpha)}{\alpha(\alpha+1)E}\eta\right]^{\frac{2\alpha-1}{\alpha-1}}~{}~{},~{}~{}0.5<\alpha<1~{},~{}\textrm{as}\
\eta\to\infty,$ (10e) and thus from (10e), we obtain $f\sim
C+\displaystyle\frac{\alpha(\alpha+1)^{2}E^{\frac{2\alpha-1}{\alpha-1}}}{-C^{2}(2\alpha-1)}\left[1+\frac{C(1-\alpha)}{\alpha(\alpha+1)E}\eta\right]^{\frac{2\alpha-1}{\alpha-1}}~{}~{},~{}~{}0.5<\alpha<1~{},~{}\textrm{as}\
\eta\to\infty.$ (10f)

By inspection, we see that (10f) approaches $C$ as $\eta\to\infty$ only when
$0.5<\alpha<1$, and thus condition (4d) can only be satisfied in this range.
Consequently, system (4d) is only valid for $0.5<\alpha\leq 1$ ($\alpha=1$ for
Newtonian fluids). Note that Fox et al. [11] incorrectly indicate that the
solution to system (4d) exists when $0<\alpha<0.5$. Additionally, although
such solutions may be obtained to the finite-domain approximation to system
(4d) where (4d) is replaced with $f^{\prime}(L)=0$, these solutions do not
converge to an infinite domain solution as $L\to\infty$; note that Pop et al.
[28] incorrectly claim that the solution exists for $\alpha>1$. The reader is
referred to comments made in Sec. I regarding the physical significance of
this restriction.

###  III.2 Construction of a Convergent Power Series Solution

To overcome the convergence limitation of the power series solution (8d), we
follow the approach of Naghshineh et al. [14] for the Newtonian problem (see
Sec. I). Here, we propose the following variable transformation, inspired by
the asymptotic expansion (10f), whose utility is validated in what follows. We
write

$U(\eta)=\left[1+\mathcal{A}\eta\right]^{\lambda},$ (11a)
$f(\eta)=F(U(\eta)),$ (11b)
$\mathcal{A}=\frac{C(1-\alpha)}{\alpha(\alpha+1)E}~{}~{},~{}~{}\lambda=\frac{2\alpha-1}{\alpha-1}.$
(11c)

The transformation (11c) maps $\eta\in[0,\infty)$ to $U\in(0,1]$ when
$0.5<\alpha<1$.

It is worth noting that the non-Newtonian transformation in (11ca), reduces to
the Newtonian transformation given by (7c), as $\alpha\to 1$, since

$\lim_{\alpha\rightarrow
1}\left[1+\mathcal{A}\eta\right]^{\lambda}=e^{\displaystyle-C\eta/2E},$ (12)

when $E=1$; we have indeed verified numerically that $E\to 1$ as $\alpha\to 1$
(see Table 2). Substituting (11c) into (4a), applying the chain rule, and
rearranging terms, we obtain the transformed ODE

$k_{1}U^{\left(3-\frac{3}{\lambda}\right)}F^{\prime\prime\prime}+k_{2}U^{\left(2-\frac{3}{\lambda}\right)}F^{\prime\prime}+k_{3}U^{\left(1-\frac{3}{\lambda}\right)}F^{\prime}-\left\\{k_{4}U^{\left(2-\frac{2}{\lambda}\right)}F^{\prime\prime}+k_{5}U^{\left(1-\frac{2}{\lambda}\right)}F^{\prime}\right\\}^{2-\alpha}F=0,$

where the primes denote derivatives of $F$ with respect to $U$. After
multiplying the above by $U^{\left(\frac{3}{\lambda}-1\right)}$, and
rearranging the ODE such that the highest derivative is on the left side of
the equation, we obtain

$k_{1}U^{2}F^{\prime\prime\prime}=-k_{2}UF^{\prime\prime}-k_{3}F^{\prime}+\left\\{k_{4}UF^{\prime\prime}+k_{5}F^{\prime}\right\\}^{2-\alpha}F,$
(13a) where $\displaystyle k_{1}=$
$\displaystyle\alpha(\alpha+1)\mathcal{A}^{3}\lambda^{3},$ (13b)
$\displaystyle k_{2}=$ $\displaystyle
3\alpha(\alpha+1)\mathcal{A}^{3}\lambda^{2}(\lambda-1),$ (13c) $\displaystyle
k_{3}=$ $\displaystyle
3\alpha(\alpha+1)\mathcal{A}^{3}\lambda(\lambda-1)(\lambda-2),$ (13d)
$\displaystyle k_{4}=$ $\displaystyle-\mathcal{A}^{2}\lambda^{2},$ (13e)
$\displaystyle k_{5}=$ $\displaystyle-\mathcal{A}^{2}\lambda(\lambda-1),$
(13f) and the expressions for $\mathcal{A}$ and $\lambda$ are defined in
(11c). The boundary conditions at $\eta=0$ from (4b) and (4c), corresponding
to $U=1$, become $F(1)=0,$ (13g) $F^{\prime}(1)=\frac{1}{\mathcal{A}\lambda},$
(13h)
$F^{\prime\prime}(1)=\frac{\kappa}{\mathcal{A}^{2}\lambda^{2}}-\frac{\lambda-1}{\mathcal{A}\lambda^{2}},$
(13i)

where the boundary condition (9) is used to obtain $F^{\prime\prime}(1)$ in
(13ii).

We next assume a solution to (13a) of the form

$F(U)=\sum_{n=0}^{\infty}{A}_{n}U^{n},~{}|U|<U_{s}(\alpha)~{},~{}E>0,$ (14)

where $U_{s}(\alpha)$ is the radius of convergence. The $E>0$ restriction
allows for the use of (10f) without introducing branch point singularities
into the $U$ domain of interest. Series (14) is readily differentiated term-
by-term to compute $F^{\prime}$, $F^{\prime\prime}$, and
$F^{\prime\prime\prime}$. After employing JCP Miller’s formula [22] and
Cauchy’s product rule [23] (see Appendixes A.1 and A.2, respectively) to re-
order the nonlinear terms in (13a), the ODE in (13a) becomes

$\sum_{n=2}^{\infty}k_{1}(n+1)n(n-1)A_{n+1}U^{n}=(-k_{3}A_{1}+A_{0}d_{0})+(-2k_{2}A_{2}-2k_{3}A_{2}+A_{0}d_{1}+A_{1}d_{0})U+\\\
\sum_{n=2}^{\infty}\left(-k_{2}(n+1)nA_{n+1}-K_{3}(n+1)+A_{0}\tilde{d}_{n}+\frac{(2-\alpha)A_{0}d_{0}}{c_{0}}(n+1)(k_{4}n+k5)A_{n+1}+\tilde{e}_{n}\right)U^{n},$
(15)

where

$\tilde{d}_{n}=\frac{1}{nc_{0}}\sum_{j=1}^{n-1}(3j-\alpha
j-n)c_{j}d_{n-j}~{}~{},~{}~{}\tilde{e}_{n}=\sum_{j=1}^{n}A_{j}d_{n-j},$ (16a)
$d_{n>0}=\frac{1}{nc_{0}}\sum_{j=1}^{n}(3j-\alpha
j-n)c_{j}d_{n-j}~{}~{},~{}~{}d_{0}=(c_{0})^{(2-\alpha)},$ (16b)
$c_{n>0}=(n+1)A_{n+1}(k_{4}n+k_{5})~{}~{},~{}~{}c_{0}=k_{5}A_{1}.$ (16c)

Using the asymptotic solution (10f), we enforce

$A_{0}=C.$ (16d)

Equating constant terms on both sides of (15) leads to

${A}_{1}=\left[\frac{k_{3}}{A_{0}(k_{5})^{2-\alpha}}\right]^{\frac{1}{1-\alpha}},$
(16e)

and equating $U^{1}$ terms on both sides of (15) leads to

$A_{2}=\frac{A_{1}d_{0}}{2k_{2}+2k_{3}-\frac{2(2-\alpha)A_{0}d_{0}}{c_{0}}(n+1)(k_{4}n+k5)}.$
(16f)

We equate like-terms in (15) to obtain the coefficients $A_{n+1}$, and for
$n\geq 2$, we obtain the recurrence relation

$A_{n+1}=\frac{A_{0}\tilde{d}_{n}+\tilde{e}_{n}}{K_{1}(n+1)n(n-1)+k_{2}(n+1)n+k_{3}(n+1)-\frac{(2-\alpha)A_{0}d_{0}}{c_{0}}(n+1)(k_{4}n+k5)}~{},~{}n\geq
2,$ (16g)

where $\tilde{d}_{n}$ and $\tilde{e}_{n}$ are defined in (16a). Transforming
back to $f(\eta)$ space via (11c), our expansion about $U=0$ (i.e.
$\eta=\infty$) is

$f(\eta)=\sum_{n=1}^{\infty}A_{n}\left[1+\mathcal{A}\eta\right]^{\lambda}.$
(16h)

Although our ultimate goal is to have a self-contained solution that is not
dependent on numerically determined parameters, at this stage we use the
numerical values of $\kappa$, $C$, and $E$ in (16h) (see Table 3) to assess
the efficacy of (16h). To determine the numerical value of $E$, we solve (10d)
for $E$ as

$E=(-h^{\prime\prime})^{\alpha-1}-\frac{C(1-\alpha)}{\alpha(\alpha+1)}\eta~{}~{}\textrm{as}\
\eta\to\infty,$ (17)

where the $\eta\to\infty$ condition is approximated in a numerical solution
(of domain length $L$) by replacing $\eta$ with $L$, $C$ with $h(L)$, and
$h^{\prime\prime}$ with $h^{\prime\prime}(L)$. Figure 2 shows that the
transformed series solution (16h) (dotted curves) matches the numerical
solution as $\eta\to\infty$, and the standard power series solution (8d)
(dashed curves) matches the numerical solution as $\eta\to 0$, as expected. It
is apparent here that the power series solution (16h) diverges as
$\eta\rightarrow 0$. Since the coefficients of (16h) alternate in sign, the
closest singularity lies along the negative real $U$ axis[29], i.e. outside
the physical domain. The vertical solid line (marked by arrow B) in Fig. 2
shows the radius of convergence of the series solution (16h). This radius is
confirmed via a numerical ratio test in the form of a Domb-Skyes plot[29],
shown in Fig. 4 as a plot of the relevant coefficient ratio vs. $1/n$. As the
curve is linear in $1/n$ for large $n$, the radius of convergence is
identified as the $y$-intercept.

Figure 4: Domb-Sykes plot for (16h) with $\alpha=0.8$, Here, the intercept for
$1/n=0$ yields the numerical radius of convergence ($y$-axis) $U_{s}\approx
0.37$, in agreement with a radius of convergence in the original domain
$\eta_{s}\approx 0.95$ through (11ca). This is consistent with the divergent
behavior observed in Fig. 2.

Naghshineh et al. [14] show that the radius of convergence for the Newtonian
Sakiadis problem can be increased by changing the expansion point of the power
series solution to the transformed ODE. Inspired from their work, we change
the expansion point in the power series solution to (13i) to $U=1$,
corresponding to $\eta=0$. Using the same procedures employed above to obtain
the series about $U=0$, the series about $U=1$ is defined as

$F(U)=\sum_{n=0}^{\infty}\hat{A}_{n}(U-1)^{n},~{}E>0,$ (18a) where
$\hat{A}_{0}=0,$ (18b) $\hat{A}_{1}=\frac{1}{\mathcal{A}\lambda},$ (18c)
$\hat{A}_{2}=\frac{1}{2}\left(\frac{\kappa}{\mathcal{A}^{2}\lambda^{2}}-\frac{\lambda-1}{\mathcal{A}\lambda^{2}}\right),$
(18d)
$\hat{A}_{3}=\frac{1}{6k_{1}}\left(-2k_{2}\hat{A}2-k_{3}\hat{A}_{1}\right),$
(18e)
$\hat{A}_{4}=\frac{1}{24k_{1}}\left(-2k_{2}\hat{A}2-6k_{2}\hat{A}_{3}-2k_{3}\hat{A}_{2}+\hat{A}_{1}\hat{d}_{0}-12k_{1}\hat{A}_{3}\right).$
(18f)

Following the same approach as was employed earlier, we obtain

$\hat{A}_{n+3}=\frac{\left\\{-k_{2}(n+1)n-k_{3}(n+1)-k_{1}(n+1)n(n-1)\right\\}\hat{A}_{n+1}}{k_{1}(n+3)(n+2)(n+1)}+\\\
\frac{\left\\{-k_{2}(n+2)(n+1)-2k_{1}(n+2)(n+1)n\right\\}\hat{A}_{n+2}+\hat{e}_{n}}{k_{1}(n+3)(n+2)(n+1)}~{},~{}n\geq
2,$ (18g)

with

$\hat{d}_{n}=\frac{1}{n\hat{c}_{0}}\sum_{j=1}^{n}(3j-\alpha
j-n)\hat{c}_{j}\hat{d}_{n-j}~{}~{},~{}~{}\hat{d}_{0}=(\hat{c}_{0})^{(2-\alpha)}\\\
$ (18h) $\hat{e}_{n}=\sum_{j=0}^{n}\hat{A}_{j}\hat{d}_{n-j}$ (18i)
$\hat{c}_{n>0}=(k_{4}n+k_{5})(n+1)\hat{A}_{n+1}+k_{4}(n+2)(n+1)\hat{A}_{n+2}~{}~{},~{}~{}\hat{c}_{0}=2k_{4}\hat{A}_{2}+k_{5}\hat{A}_{1}.$
(18j)

The constants $k_{1}$ \- $k_{5}$ are defined in (13b) - (13f).

Figure 2 shows that the transformed series solution (18j) (solid curves)
matches the numerical solution as $\eta\to\infty$, as well as $\eta\to 0$.
Thus, moving the location of the expansion point from $U=0$ to $U=1$ enables a
convergent expansion over the whole domain. It should be noted here that this
is distinctly different from the Newtonian case[14], where both expansions
about $U=0$ and $U=1$ converge–although the latter expansion converges faster.
Figure 5a shows the absolute error (the absolute difference) between $N$-term
truncations of the convergent series solution (18j) and the numerical solution
for $\alpha=0.8$. The numerical values of the constants used to generate this
figure are shown in Appendix C.2. Here, we choose to stop at $N=200$ in
construction of Fig. 5 because the absolute error is close to machine
precision. The dashed curve in Fig. 6 shows the infinity norm (maximum
absolute error) between the $N$-term truncation of the series solution (18j)
and the numerical solution taken over $\eta\in[0,L]$ with $L=11000$ for
$\alpha=0.8$, using the values of constants generated by the numerical
solution (see Appendix C.2 for more details). The plateau reached in this
figure occurs when further refinements to the power series lead to errors
smaller than that of the numerical solution.

(a)  (b)

Figure 5: Absolute error (the absolute value of the difference) between
$N$-term truncations of the convergent series solution (18j) and the numerical
solution (over the domain $\eta\in[0,L]$ with $L=11000$) for $\alpha=0.8$,
plotted versus $\eta$, using $\kappa$, $C$, and $E$ values (a) generated by
the numerical solution, and (b) predicted by equations system (19c), discussed
in Sec. III.3.  (a)  Figure 6: The infinity norm (maximum absolute error)
between $N$-terms truncations of (18j) and the numerical solution (occurring
over $\eta\in[0,L]$ with $L=11000$) for $\alpha=0.8$, plotted versus $N$. The
dashed and solid curves correspond to when numerical and algorithmically
predicted values (using equations system (19c) discussed in Sec. III.3) of
$\kappa$, $C$, and $E$ are used, respectively. Convergence of the numerical
and algorithmically predicted values of the constants are reported in
respective Tables 3 and 4 of Appendix C.2.

###  III.3 Prediction of Unknown Parameters

Thus far, results have been presented where the numerically-obtained values of
the constants $\kappa$, $C$, and $E$ have been used. With the aim of making
the series solution (18j) independent of the numerical solution, we adapt an
algorithm used by Barlow et al. [13] to predict the values of the constants
$\kappa$, $C$, and $E$. In that study, it is sufficient to construct a system
of equations for the unknowns by choosing the last $\mathcal{N}$ coefficients
of the series solution to be zero, where $\mathcal{N}$ is the number of
unknowns ($\mathcal{N}=3$ in this problem: $\kappa$, $C$, and $E$). For a
convergent series, it is most certainly the case that the last coefficient
approaches zero as $N\to\infty$, so this assumption is self consistent in the
limit. Additionally, from the perspective of the number of equations and
unknowns, we need three equations to find $\kappa$, $C$, and $E$. The
assumption of series convergence implicit in the equations to solve is
validated by convergence of the algorithm itself for increasingly large
numbers of series terms. However, in this non-Newtonian extension we find that
the three equations ($\hat{A}_{N}=0$, $\hat{A}_{N-1}=0$, and
$\hat{A}_{N-2}=0$) are linearly dependent, as evident by the determinant of
the Jacobian. For that reason, we alter the algorithm such that we use one of
those equations ($\hat{A}_{N}=0$) (19c). For the remaining two equations we
adapt the algorithm used by Naghshineh et al. [14] where conditions are
imposed from the side of the domain that is opposite to that of the series’
expansion point. The system of equations used here are:

$\left[\sum_{n=0}^{N}\hat{A}_{n}(-1)^{n}\right]-C=0,$ (19a)
$\left[\sum_{n=0}^{N}n\hat{A}_{n}(-1)^{n-1}\right]-A_{1}=0,$ (19b)
$\hat{A}_{N}=0,$ (19c)

where $C$, $A_{1}$, and $\hat{A}_{N}$ are defined in (6), (16e), and (18jg),
respectively. Note that (19a) and (19b) correspond to the boundary conditions
$F(U=0)=C$ and $F^{\prime}(U=0)=A_{1}$. Newton’s Method is used to solve the
system of equations (19c), and the details are provided in Appendix D.

Figures 5b and 6 (solid curve) show typical results of the solution of
equations system (19c), with a tolerance of $10^{-15}$ used in the
implementation of Newton’s method. Figure 5b shows the absolute error vs
$\eta$ compared with the numerical solution when $\kappa$, $C$, and $E$ are
predicted algorithmically for $\alpha=0.8$. The solid curve in Fig. 6 shows
the maximum absolute error vs $N$ compared with the numerical solution when
$\kappa$, $C$, and $E$ are predicted using the solution of equations system
(19c). As seen by inspection, the accuracy of the solution, $f(\eta)$,
increases with the number of terms used in the series. It is important to note
that we do not explicitly enforce that $f(\eta\rightarrow\infty)\to C$ (for
all $N$) at $U=0$ (corresponding to $\eta\rightarrow\infty$) in (18ja), as we
do for the expansion in (14) (refer to (16d)). Consequently, for $N=100$ and
200, this allows for the curves in Fig. 5b to ultimately attain lower error
values than those shown in Fig. 5a.

The key issue with using the series solution (18j) is its need to access the
value of the asymptotic constant $E$. Although not shown here, we have
examined a variety of other permissible $\alpha$ values, and we find that the
numerical $E$ ultimately becomes negative for $\alpha<0.74$, thus invalidating
the use of the gauge function (11c) due to a branch point singularity that
arises in the physical $\eta$ domain; in this range of $\alpha$, equations
system (19c) fails to predict converged values for $\kappa$, $C$, and $E$.
That said, the asymptotic form with $E$ being negative is perfectly valid for
large enough $\eta$, which is the region in which the asymptotic form itself
is valid.

##  IV Asymptotically Motivated Approximant

Since a convergent power series solution has only been obtained for
$0.74<\alpha<1$, we consider an alternative approach to obtain an analytical
form over the full range of $0.5<\alpha<1$. Note that the form we will obtain
can be used as an alternative to the convergent series over the full range of
$\alpha$, although some numerical results are needed to do so. One way of
overcoming convergence barriers in divergent series solutions is to
analytically continue them via Padé approximants [30]. To this end, we utilize
an asymptotically motivated approximant[13] in the form of a modified Padé
approximant as

$f_{A}=C-\left[\frac{\displaystyle\sum_{n=0}^{M+1}P_{n}\eta^{n}}{\displaystyle\sum_{n=0}^{M}Q_{n}\eta^{n}}\right]^{\lambda}.$
(20a) In (20ba), $C$ and $\lambda$ are given by (6) and (11c), respectively.
To solve for the coefficients $P_{n}$ and $Q_{n}$ in (20ba), we write (20ba)
in the form
$\left[C-\sum_{n=0}^{\infty}a_{n}\eta^{n}\right]^{1/\lambda}=~{}\frac{\displaystyle\sum_{n=0}^{M+1}P_{n}\eta^{n}}{\displaystyle\sum_{n=0}^{M}Q_{n}\eta^{n}},$
(20b)

where the coefficients $a_{n}$ are given in (8da). Using JCP Miller’s formula
(see Appendix A.1) on the left hand side of (20bb), a standard Padé solver may
be employed to solve for the coefficients $P_{n}$ and $Q_{n}$ on the right-
hand side of (20bb). Note that for $\eta\rightarrow\infty$, $f_{A}\sim
C-(P_{M+1}/Q_{M})^{\lambda}\eta^{\lambda}$, which is consistent with the
asymptotic form (10f) as $\eta\gg E$. To implement the approximant, we use the
numerically predicted values of $\kappa$ and $C$, reported in Table 1. This is
precisely the approach taken by Belden et al. [21] in the asymptotic
approximant solution to the Falkner Skan equation. More details about the
numerical prediction for various domain lengths $L$ for $\alpha=0.8$ are
provided in Appendix C.2. As seen in Table 1, the numerical solution of the
boundary value problem reveals a high sensitivity of parameter values to
domain length, and this sensitivity increases as $\alpha$ decreases. Figures 7
and 8 show the absolute error between $M$-term truncations of approximant
(20b) and the numerical solution for $\alpha=0.8$ and $\alpha=0.6$,
respectively.

Table 1: Numerical values of the constants $\kappa$ and $C$ for the non-Newtonian Sakiadis problem. The values are computed using the shooting method algorithm explained in Appendix B with $\eta=\infty$ replaced with a finite surrogate $L$ given in the table. The values of $\kappa$ and $C$ are accurate to within the decimal places reported here, based on convergence by successively increasing $L$. $\alpha$ | $L$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\kappa$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}C$
---|---|---|---
0.99 | 100 | -0.4434518189261262 | 1.6250769221853265
0.9 | 2000 | -0.4413601253597191 | 1.717915813011322
0.8 | 11000 | -0.440672715940425 | 1.860152537
0.7 | 11000 | -0.442664523 | 2.08739
0.6 | 11000 | -0.44906693 | 2.56
0.55 |  40000 | -0.454851 | 3.1

There are a few defective approximants 333It is possible that poles of a Padé
approximant arise within the physical domain of a problem for particular
degrees of denominator and numerator. If the exact solution is expected to be
finite within the physical domain, these Padés are deemed defective. For this
exact reason, we did not use $M=25$ in Fig. 8 as it led to a defective
approximant for the case of $\alpha=0.6$. that arise between the indicated
truncations of the approximant in Figs. 7 and 8, in which the denominator in
(20bb) becomes zero for positive $\eta$ values; it is standard practice to
ignore these when assessing the solution [30]. We note that the smallest error
is obtained at $M=26$ and $M=15$ for the cases of $\alpha=0.8$ and $0.6$,
respectively. For any larger value of $M$, the error oscillates between curves
that are similar to $M=25$ and $27$ in Fig. 7 and $M=12$ and $26$ in Fig. 8.
In both cases, we accept the solution as converged, as its precision (defined
here as amplitude of these oscillations) is consistent with that of the
inputs, particularly $C$ (see Table 1). Figure 9 shows the infinity norm
(maximum absolute error) vs $\alpha$ for the permissible range of $\alpha$
values when $M=20$, using the numerical values of the constants $\kappa$ and
$C$.

(a)  Figure 7: Absolute error between $M$-term truncations of the approximant
(20b) and the numerical solution (occurring over $\eta\in[0,L]$ with
$L=11000$) for $\alpha=0.8$, plotted versus $\eta$, using $\kappa$ and $C$
generated by the numerical solution. Error decreases until $M=26$, after which
it oscillates between curves that are similar to $M=25$ and 27. (a)  Figure 8:
Absolute error between $M$-terms truncations of the approximant (20b) and the
numerical solution (occurring over $\eta\in[0,L]$ with $L=11000$) for
$\alpha=0.6$, plotted versus $\eta$, using $\kappa$ and $C$ generated by the
numerical solution. Error decreases until $M=15$, after which it oscillates
between curves that are similar to $M=12$ and 26. (a)  Figure 9: Maximum
absolute error between approximant (20b) (with $M=20$) and the numerical
solution (occurring over $\eta\in[0,L]$) plotted versus $\alpha$, using
$\kappa$ and $C$ generated by the numerical solution. The values of $C$ and
$\kappa$ are shown in Table 1.

##  V Post-Processing: Analytically obtained streamlines

Now that we have accurate analytical solutions to (4d), we may insert $f$,
given by either (16h) (for $0.74\leq\alpha\leq 1$)444For $\alpha=1$, one can
use the Newtonian result (7c) in place of (16h); coefficients for (7c) are
provided by Naghshineh et al[14]. or (20b) (for $0.5<\alpha<1$), and its
derivative $f^{\prime}(\eta)$ (which may obtained analytically) into (3) to
obtain the $u$ velocity field, which is shown in the right-hand plot in Fig.
10. The paths of fluid points, i.e., the streamlines of constant $\psi$, may
be extracted easily from the analytical solution. To do so, we explicitly
solve for the $x$ and $y$ coordinates of a given streamline $\psi=$ constant
by rearranging the equations in (2b) to yield:

$y=\frac{\psi\eta}{u_{w}f(\eta)},~{}~{}x=\frac{\rho}{Ku_{w}^{2\alpha-1}}\left(\frac{\psi}{f(\eta)}\right)^{\alpha+1}.$
(21)

Equation (21) provides a parametric representation of the streamlines in terms
of $\eta$ and $f(\eta)$, the latter given analytically by (16h) or (20b).
Figure 10 provides a typical streamline plot extracted in this way. In the
figure, a dashed curve plots the boundary layer thickness
$y=\eta\left(K~{}x~{}u_{w}^{\alpha-2}/\rho\right)^{1/(\alpha+1}$, defined here
as the locus of points where the fluid velocity is reduced to 10% of the wall
velocity; according to (2a), this occurs when $u/u_{w}=df/d\eta=0.1$. From the
topmost plot of Fig. 10, this occurs when $\eta=4.04$. The benefit of the
analytical solution is clearly indicated here, as streamline plots can be
generated accurately to any desired resolution with low computational cost.

(a)  Figure 10:  (top) $df/d\eta$ obtained analytically from (16h) for
$\alpha=0.8$ (using $C$, $\kappa$, and $E$ values from Table 4 with $N=400$,
using system of equations (19c)), with gridlines indicating that
$df/d\eta=u/u_{w}=0.1$ at $\eta\approx 4.04$. (left) Contours of constant
$\psi$ obtained analytically from (16h) (using 27 terms) and (2b), displayed
in increments of $\Delta\psi=1$ in the $\bar{y}\equiv u_{w}y$ vs.
$\bar{x}\equiv Ku_{w}^{2\alpha-1}x/\rho$ plane. The dashed curve is the
envelope of the boundary layer (chosen here to be the locus of points at which
$u/u_{w}=0.1$), and restricts the display of the streamlines for velocities
where $u/u_{w}>0.1$. (right) The $u$ velocity field obtained analytically from
(16h) (using 27 terms) and (3).

##  VI Conclusions

In this work, we provide a convergent power series solution to the non-
Newtonian Sakiadis boundary layer problem, valid for $0.74\leq\alpha<1$, using
the asymptotic expansion as $\eta\to\infty$ to determine a gauge function for
the series. The asymptotically motivated series fails when the gauge function
is unable to completely transverse the physical domain due to a branch point
singularity that arises in the asymptotic form. We note that although we
developed an asymptotic approximant to model cases where $0.5<\alpha<0.74$,
the approximant is capable of representing the solution for all $\alpha$
values in the range $0.5<\alpha<1$. Once obtained, the analytical solutions
enable computationally-efficient post-processing to extract streamlines to any
desired resolution.

##  Appendix A Useful Formulae for Manipulating Series

###  A.1 Raising a series to a power

The following relation is JCP Miller’s formula for raising a series to a power
[22]:

$\left(\sum_{n=0}^{\infty}a_{n}x^{n}\right)^{\gamma}=\sum_{n=0}^{\infty}b_{n}x^{n},$
(22a) where
$b_{n>0}=\frac{1}{n~{}a_{0}}\sum_{j=1}^{n}(j\gamma-n+j)a_{j}b_{n-j}~{},~{}b_{0}=\left(a_{0}\right)^{\gamma}~{},~{}a_{0}\neq
0.$ (22b)

###  A.2 Product of two series

The following relation is the well-known Cauchy product of two series [23]:

$\sum_{n=0}^{\infty}a_{n}x^{n}\sum_{n=0}^{\infty}b_{n}x^{n}=\sum_{n=0}^{\infty}\left(\sum_{j=0}^{n}a_{j}b_{n-j}\right)x^{n}.$
(23)

##  Appendix B  Numerical Solution: Shooting Method

The algorithm below is motivated from the work Cebeci et al. [20], developed
for the Falkner-Skan boundary layer problem. Here, we extend their approach to
the non-Newtonian Sakiadis problem given by (4d). With the goal of determining
the value of $\kappa$ defined by (9), we approximate the boundary value
problem (4d) in $f(\eta)$ (defined on a semi-infinite domain) with the
following initial value problem (IVP) in $f(\eta;\kappa)$ (defined on a finite
domain):

$\alpha(\alpha+1)f^{\prime\prime\prime}-(-f^{\prime\prime})^{(2-\alpha)}f=0,~{}~{}0\leq\eta\leq
L$ (24a) with the initial conditions (taken from (4b), (4c), and (9))
$f(0;\kappa)=0~{}~{},~{}~{}f^{\prime}(0;\kappa)=1~{}~{},~{}~{}f^{\prime\prime}(0;\kappa)=\kappa,$
(24b)

where $f^{\prime}$ denotes the derivative of $f$ with respect to $\eta$. In
order to determine $\kappa$, we subject the IVP (24b) to the constraint

$f^{\prime}(L;\kappa)=0,$ (25)

which incorporates condition (4d) such that the solution to (24b) limits to
the solution of (4d) as $L\to\infty$. The determination of $\kappa$ and the
numerical solution for $f$ itself from (24b) (subject to (25)) is obtained by
the method of shooting, as outlined below.

First, we replace (24a) with a system of three first-order ODEs. To do so, we
let $\mathcal{U}(\eta;\kappa)$, $\mathcal{V}(\eta;\kappa)$, and
$\mathcal{V}^{\prime}(\eta;\kappa)$ represent $f^{\prime}(\eta;\kappa)$,
$f^{\prime\prime}(\eta;\kappa)$, and $f^{\prime\prime\prime}(\eta;\kappa)$,
respectively. Thus, (24b) can be written as the system

$f^{\prime}=\mathcal{U}~{}~{},~{}~{}\mathcal{U}^{\prime}=\mathcal{V}~{}~{},~{}~{}\mathcal{V}^{\prime}=\frac{(-\mathcal{V})^{(2-\alpha)}f}{\alpha(\alpha+1)},$
(26a) with initial conditions
$f(0;\kappa)=0~{}~{},~{}~{}\mathcal{U}(0;\kappa)=1~{}~{},~{}~{}\mathcal{V}(0;\kappa)=\kappa.$
(26b)

The objective is to provide a solution to the IVP (26), such that the
constraint (25) is satisfied. That is, we solve for the solution of (26) by
seeking $\kappa$, such that

$\mathcal{U}(L;\kappa)=0,$ (27)

where $L$ is successively increased, such that the $\kappa$ value for (24b)
approaches the $\kappa$ value for (4d).

In order to determine $\kappa$ in (26b), we use Newton’s method [33] defined
by

$\kappa^{\gamma+1}=\kappa^{\gamma}-\frac{\mathcal{U}(L;\kappa^{\gamma})}{\frac{\partial}{\partial\kappa}(\mathcal{U}(L;\kappa^{\gamma}))}~{}~{},~{}~{}\gamma=0,~{}1,~{}2,~{}...,$
(28)

with $\kappa^{0}$ being the initial estimate for $\kappa$, and $\gamma$ is the
iteration number. In order to obtain the derivative of $\mathcal{U}$ with
respect to $\kappa$, we take the derivative of (26) with respect to $\kappa$,
which leads to the following additional IVP:

$\frac{\partial
f^{\prime}}{\partial\kappa}=\frac{\partial\mathcal{U}}{\partial\kappa}~{}~{},~{}~{}\frac{\partial\mathcal{U}^{\prime}}{\partial\kappa}=\frac{\partial\mathcal{V}}{\partial\kappa}~{}~{},~{}~{}\frac{\partial\mathcal{V}^{\prime}}{\partial\kappa}=(\alpha-2)\frac{\partial\mathcal{V}}{\partial\kappa}~{}\frac{(-\mathcal{V})^{(1-\alpha)}}{\alpha(\alpha+1)}f+\frac{(-\mathcal{V})^{(2-\alpha)}}{\alpha(\alpha+1)}~{}\frac{\partial
f}{\partial\kappa},$ (29a) with conditions $\frac{\partial
f}{\partial\kappa}(0;\kappa)=0~{}~{},~{}~{}\frac{\partial\mathcal{U}}{\partial\kappa}(0;\kappa)=0~{}~{},~{}~{}\frac{\partial\mathcal{V}}{\partial\kappa}(0;\kappa)=1.$
(29b)

For clarity, we assign new variables to the derivatives with respect to
$\kappa$ as follows

$\mathfrak{F}(\eta;\kappa)\equiv\frac{\partial
f}{\partial\kappa}(\eta;\kappa),~{}~{}\mathfrak{U}(\eta;\kappa)\equiv\frac{\partial\mathcal{U}}{\partial\kappa}(\eta;\kappa),~{}~{}\mathfrak{V}(\eta;\kappa)\equiv\frac{\partial\mathcal{V}}{\partial\kappa}(\eta;\kappa),$

commute the differentiation with respect to $\eta$ (denoted by primes) and
differentiation with respect to $\kappa$ in (29b), and combine (26) and (29b)
into the single IVP evaluated at $\kappa=\kappa^{\gamma}$:

$\frac{d}{d\eta}\left[\begin{array}[]{ccccccc}f\\\ \mathcal{U}\\\
\mathcal{V}\\\ \mathfrak{F}\\\ \mathfrak{U}\\\
\mathfrak{V}\end{array}\right]=\left[\begin{array}[]{ccccccc}\mathcal{U}\\\
\mathcal{V}\\\ \frac{(-\mathcal{V})^{(2-\alpha)}f}{\alpha(\alpha+1)}\\\
\mathfrak{U}\\\ \mathfrak{V}\\\
(\alpha-2)\mathfrak{V}\frac{(-\mathcal{V})^{(1-\alpha)}}{\alpha(\alpha+1)}f+\frac{(-\mathcal{V})^{(2-\alpha)}}{\alpha(\alpha+1)}\mathfrak{F}\end{array}\right],~{}~{}\left[\begin{array}[]{ccccccc}f\\\
\mathcal{U}\\\ \mathcal{V}\\\ \mathfrak{F}\\\ \mathfrak{U}\\\
\mathfrak{V}\end{array}\right]_{\eta=0}=\left[\begin{array}[]{ccccccc}0\\\
1\\\ \kappa^{\gamma}\\\ 0\\\ 0\\\ 1\end{array}\right].$ (30)

From (30), both $\mathcal{U}(L;\kappa^{\gamma})$ and
$\mathfrak{U}(L;\kappa^{\gamma})$ may be determined; then the guess
$\kappa^{\gamma}$ may be progressed to the next iteration via

$\kappa^{\gamma+1}=\kappa^{\gamma}-\frac{\mathcal{U}(L;\kappa^{\gamma})}{\mathfrak{U}(L;\kappa^{\gamma})}~{}~{},~{}~{}\gamma=0,~{}1,~{}2,~{}...,$
(31)

to compute $\kappa^{\gamma+1}$. We use a fourth-order Runge-Kutta method, with
$\Delta\eta=0.001$ to solve (30). A convergence requirement of
$|\kappa^{\gamma+1}-\kappa^{\gamma}|<10^{-15}$ is enforced in the Newton
iteration (31).

##  Appendix C  Evaluations of $\kappa$, $C$, and $E$

###  C.1 Numerical evaluation of $E$

In this section, we show that the constant $E\to 1$ as $\alpha\to 1$. These
values are obtained numerically via a shooting method described in Appendix B
combined with the definition of $E$ defined in (17).

Table 2: Numerical values of the constant $E$ defined in (17). The values are computed using the shooting method algorithm explained in Appendix B with $\eta=\infty$ replaced with a finite surrogate $L$ given in the table.  $\alpha$ | $L$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}E$
---|---|---
0.99 |  100 | 0.996357144778925
0.9 |  2000 | 0.921210805240775
0.8 |  11000 | 0.639133894794667
0.75 |  11000 | 0.249305949228073

###  C.2  Evaluations of $\kappa$, $C$, and $E$ for $\alpha=0.8$

In what follows, we provide constants used in the analytical solutions (8d),
(11c), (18j), and (20b) for $\alpha=0.8$. These are obtained in two
ways–either numerically via a shooting method (Table 3 and Appendix B), or via
a self-contained algorithm using the power series (Table 4 and Appendix (D)).

Table 3: Numerical values of the constants $\kappa$, $C$, and $E$ with $\alpha=0.8$ for the non-Newtonian Sakiadis problem. The values are computed using the shooting method algorithm explained in Appendix B with $\eta=\infty$ replaced with a finite surrogate $L$ given in the table. This numerical algorithm uses a Newton iteration with a tolerance of $10^{-15}$, wrapped around a 4th-order Runge-Kutta solver with a step size of $\triangle\eta=0.001$.  $L$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\kappa$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}C$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}E$
---|---|---|---
40 | -0.44069277587148 | 1.85680030949222 | 0.654388901675618
100 | -0.44067330624905 | 1.85989634054125 | 0.641907516688146
240 | -0.44067273476661 | 1.86013263767144 | 0.639639542413533
540 | -0.44067271669187 | 1.86015073982904 | 0.639235534657387
1000 | -0.44067271600486 | 1.86015225132246 | 0.639163575176951
2000 | -0.44067271594447 | 1.86015250138706 | 0.639141167249476
5000 | -0.44067271594052 | 1.86015253504567 | 0.639134860809690
11000 | -0.44067271594043 | 1.86015253716140 | 0.639133894794667
Table 4: Predicted values of the constants $\kappa$, $C$, and $E$ with $\alpha=0.8$ for the non-Newtonian Sakiadis problem. The values are computed using the the solution of equations system (19c) explained in Sec.(III.3) with the $N$ values given in the table. This algorithm uses a Newton iteration with a tolerance of $10^{-15}$. $N$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\kappa$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}C$ | $~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}E$
---|---|---|---
50 | -0.440672770686287 | 1.860151847487335 | 0.639137961392650
100 | -0.440672715934279 | 1.860152537361902 | 0.639133656731447
200 | -0.440672715934271 | 1.860152537362062 | 0.639133656729539
400 | -0.440672715934270 | 1.860152537362068 | 0.639133656729520

##  Appendix D  Newton’s Method Used to Predict Constants

The system of equations used in this algorithm are:

$\mathcal{F}_{1}(\kappa,C,E)=\left[\sum_{n=0}^{N}\hat{A}_{n}(-1)^{n}\right]-C,$
(32a)
$\mathcal{F}_{2}(\kappa,C,E)=\left[\sum_{n=0}^{N}n\hat{A}_{n}(-1)^{n-1}\right]-A_{1},$
(32b) $\mathcal{F}_{3}(\kappa,C,E)=\hat{A}_{N},$ (32c)

where $C$, $A_{1}$, and $\hat{A}_{N}$ are defined in (6), (16e), and (18jg),
respectively.

We take the derivatives of (32a)-(32c) with respect to $\kappa$, $C$, and $E$,
and obtain the Jacobian matrix,

$\textbf{J}=\left[\begin{array}[]{ccccccc}\frac{\partial\mathcal{F}_{1}}{\partial
C}&&&\frac{\partial\mathcal{F}_{1}}{\partial
E}&&&\frac{\partial\mathcal{F}_{1}}{\partial\kappa}\\\
\frac{\partial\mathcal{F}_{2}}{\partial
C}&&&\frac{\partial\mathcal{F}_{2}}{\partial
E}&&&\frac{\partial\mathcal{F}_{2}}{\partial\kappa}\\\
\frac{\partial\mathcal{F}_{3}}{\partial
C}&&&\frac{\partial\mathcal{F}_{3}}{\partial
E}&&&\frac{\partial\mathcal{F}_{3}}{\partial\kappa}\end{array}\right],$ (33)

which is used in Newton’s Method

$\left[\begin{array}[]{ccccccc}C^{j+1}\\\ E^{j+1}\\\
\kappa^{j+1}\end{array}\right]=\left[\begin{array}[]{ccccccc}C^{j}\\\ E^{j}\\\
\kappa^{j}\end{array}\right]-\textbf{J}^{-1}\left[\begin{array}[]{ccccccc}\mathcal{F}_{1}(C^{j},E^{j},\kappa^{j})\\\
\mathcal{F}_{2}(C^{j},E^{j},\kappa^{j})\\\
\mathcal{F}_{3}(C^{j},E^{j},\kappa^{j})\\\ \end{array}\right].$ (34)

To construct the Jacobian matrix (33), we take the derivative of (32a \-
(32c), with respect to $C$, $E$, and $\kappa$, as follows:

$\frac{\partial\mathcal{F}_{1}}{\partial
C}=\left[\sum_{n=0}^{\infty}\frac{\partial\hat{A}_{n}}{\partial
C}(-1)^{n}\right]-1~{}~{},~{}~{}\frac{\partial\mathcal{F}_{1}}{\partial
E}=\sum_{n=0}^{\infty}\frac{\partial\hat{A}_{n}}{\partial
E}(-1)^{n}~{}~{},~{}~{}\frac{\partial\mathcal{F}_{1}}{\partial\kappa}=\sum_{n=0}^{\infty}\frac{\partial\hat{A}_{n}}{\partial\kappa}(-1)^{n},$
(35a) $\frac{\partial\mathcal{F}_{2}}{\partial
C}=\left[\sum_{n=0}^{\infty}n\frac{\partial\hat{A}_{n}}{\partial
C}(-1)^{n-1}\right]-\frac{\partial A_{1}}{\partial C},$ (35b)
$\frac{\partial\mathcal{F}_{2}}{\partial
E}=\left[\sum_{n=0}^{\infty}n\frac{\partial\hat{A}_{n}}{\partial
E}(-1)^{n-1}\right]-\frac{\partial A_{1}}{\partial E},$ (35c)
$\frac{\partial\mathcal{F}_{2}}{\partial\kappa}=\left[\sum_{n=0}^{\infty}n\frac{\partial\hat{A}_{n}}{\partial\kappa}(-1)^{n-1}\right]-\frac{\partial
A_{1}}{\partial\kappa},$ (35d) $\frac{\partial\mathcal{F}_{3}}{\partial
C}=\frac{\partial\hat{A}_{N}}{\partial
C}~{}~{},~{}~{}\frac{\partial\mathcal{F}_{3}}{\partial
E}=\frac{\partial\hat{A}_{N}}{\partial
E}~{}~{},~{}~{}\frac{\partial\mathcal{F}_{3}}{\partial\kappa}=\frac{\partial\hat{A}_{N}}{\partial\kappa}.$
(35e)

We take the derivative of (11c) with respect to the unknowns $C$, $E$, and
$\kappa$. It is important to note that $\lambda$ is not function of any these
three constants. Hence, we only take the derivative of $\mathcal{A}$ in (11c)
as follows:

$\frac{\partial\mathcal{A}}{\partial
C}=\frac{1-\alpha}{\alpha(\alpha+1)E}~{}~{},~{}~{}\frac{\partial\mathcal{A}}{\partial
E}=\frac{-C(1-\alpha)}{\alpha(\alpha+1)E^{2}}~{}~{},~{}~{}\frac{\partial\mathcal{A}}{\partial\kappa}=0.$
(36)

We use the chain rule to take the derivatives of (13b)-(13f), as follows:

$\frac{\partial k_{1}}{\partial
C}=3\alpha(\alpha+1)\lambda^{3}\mathcal{A}^{2}\frac{\partial\mathcal{A}}{\partial
C}~{}~{},~{}~{}\frac{\partial k_{1}}{\partial
E}=3\alpha(\alpha+1)\lambda^{3}\mathcal{A}^{2}\frac{\partial\mathcal{A}}{\partial
E}~{}~{},~{}~{}\frac{\partial k_{1}}{\partial\kappa}=0,$ (37a) $\frac{\partial
k_{2}}{\partial
C}=9\alpha(\alpha+1)\lambda^{2}(\lambda-1)\mathcal{A}^{2}\frac{\partial\mathcal{A}}{\partial
C}~{}~{},~{}~{}\frac{\partial k_{2}}{\partial
E}=9\alpha(\alpha+1)\lambda^{2}(\lambda-1)\mathcal{A}^{2}\frac{\partial\mathcal{A}}{\partial
E}~{}~{},~{}~{}\frac{\partial k_{2}}{\partial\kappa}=0,$ (37b) $\frac{\partial
k_{3}}{\partial
C}=3\alpha(\alpha+1)\lambda(\lambda-1)(\lambda-2)\mathcal{A}^{2}\frac{\partial\mathcal{A}}{\partial
C},$ (37c) $\frac{\partial k_{3}}{\partial
E}=3\alpha(\alpha+1)\lambda(\lambda-1)(\lambda-2)\mathcal{A}^{2}\frac{\partial\mathcal{A}}{\partial
E}~{}~{},~{}~{}\frac{\partial k_{3}}{\partial\kappa}=0,$ (37d) $\frac{\partial
k_{4}}{\partial C}=-2\lambda^{2}\mathcal{A}\frac{\partial\mathcal{A}}{\partial
C}~{}~{},~{}~{}\frac{\partial k_{4}}{\partial
E}=-2\lambda^{2}\mathcal{A}\frac{\partial\mathcal{A}}{\partial
E}~{}~{},~{}~{}\frac{\partial k_{4}}{\partial\kappa}=0,$ (37e) $\frac{\partial
k_{5}}{\partial
C}=-2\lambda(\lambda-1)\mathcal{A}\frac{\partial\mathcal{A}}{\partial
C}~{}~{},~{}~{}\frac{\partial k_{5}}{\partial
E}=-2\lambda(\lambda-1)\mathcal{A}\frac{\partial\mathcal{A}}{\partial
E}~{}~{},~{}~{}\frac{\partial k_{5}}{\partial\kappa}=0.$ (37f)

We take derivatives of $\hat{A}_{0}$ through $\hat{A}_{4}$ in (18jb)-(18jf),
as follows

$\frac{\partial\hat{A}_{0}}{\partial
C}=0~{}~{},~{}~{}\frac{\partial\hat{A}_{0}}{\partial
E}=0~{}~{},~{}~{}\frac{\partial\hat{A}_{0}}{\partial\kappa}=0.$ (38a)
$\frac{\partial\hat{A}_{1}}{\partial
C}=\frac{-1}{\lambda\mathcal{A}^{2}}\frac{\partial\mathcal{A}}{\partial
C}~{}~{},~{}~{}\frac{\partial\hat{A}_{1}}{\partial
E}=\frac{-1}{\lambda\mathcal{A}^{2}}\frac{\partial\mathcal{A}}{\partial
E}~{}~{},~{}~{}\frac{\partial\hat{A}_{1}}{\partial\kappa}=0,$ (38b)
$\frac{\partial\hat{A}_{2}}{\partial
C}=\frac{1}{2\lambda^{2}}\left(\frac{-2\kappa}{a^{3}}+\frac{\lambda-1}{\mathcal{A}^{2}}\right)\frac{\partial\mathcal{A}}{\partial
C},$ (38c) $\frac{\partial\hat{A}_{2}}{\partial
E}=\frac{1}{2\lambda^{2}}\left(\frac{-2\kappa}{a^{3}}+\frac{\lambda-1}{\mathcal{A}^{2}}\right)\frac{\partial\mathcal{A}}{\partial
E}~{}~{},~{}~{}\frac{\partial\hat{A}_{2}}{\partial\kappa}=\frac{1}{2\lambda^{2}\mathcal{A}^{2}}.$
(38d) $\frac{\partial\hat{A}_{3}}{\partial C}=u_{1}\frac{\partial
v_{1}}{\partial C}+v_{1}\frac{\partial u_{1}}{\partial
C}~{}~{},~{}~{}\frac{\partial\hat{A}_{3}}{\partial E}=u_{1}\frac{\partial
v_{1}}{\partial E}+v_{1}\frac{\partial u_{1}}{\partial
E}~{}~{},~{}~{}\frac{\partial\hat{A}_{3}}{\partial\kappa}=u_{1}\frac{\partial
v_{1}}{\partial\kappa},$ (38e) where
$u_{1}=\frac{1}{6k_{1}}~{}~{},~{}~{}v_{1}=-2k_{2}\hat{A}_{2}-k_{3}\hat{A}_{1},$
(38f) $\frac{\partial u_{1}}{\partial C}=\frac{-1}{6(k_{1})^{2}}\frac{\partial
k_{1}}{\partial C}~{}~{},~{}~{}\frac{\partial u_{1}}{\partial
E}=\frac{-1}{6(k_{1})^{2}}\frac{\partial k_{1}}{\partial
E}~{}~{},~{}~{}\frac{\partial u_{1}}{\partial\kappa}=0,$ (38g) $\frac{\partial
v_{1}}{\partial C}=-2\left(k_{2}\frac{\partial\hat{A}_{2}}{\partial
C}+\hat{A}_{2}\frac{\partial k_{2}}{\partial
C}\right)-\left(k_{3}\frac{\partial\hat{A}_{1}}{\partial
C}+\hat{A}_{1}\frac{\partial k_{3}}{\partial C}\right),$ (38h) $\frac{\partial
v_{1}}{\partial E}=-2\left(k_{2}\frac{\partial\hat{A}_{2}}{\partial
E}+\hat{A}_{2}\frac{\partial k_{2}}{\partial
E}\right)-\left(k_{3}\frac{\partial\hat{A}_{1}}{\partial
E}+\hat{A}_{1}\frac{\partial k_{3}}{\partial E}\right),$ (38i) $\frac{\partial
v_{1}}{\partial\kappa}=-2k_{2}\frac{\partial\hat{A}_{2}}{\partial\kappa}.$
(38j) $\frac{\partial\hat{A}_{4}}{\partial C}=u_{2}\frac{\partial
v_{2}}{\partial C}+v_{2}\frac{\partial u_{2}}{\partial
C}~{}~{},~{}~{}\frac{\partial\hat{A}_{4}}{\partial E}=u_{2}\frac{\partial
v_{2}}{\partial E}+v_{2}\frac{\partial u_{2}}{\partial
E}~{}~{},~{}~{}\frac{\partial\hat{A}_{4}}{\partial\kappa}=u_{2}\frac{\partial
v_{2}}{\partial\kappa},$ (38k) and
$u_{2}=\frac{1}{24k_{1}}~{}~{},~{}~{}v_{2}=-2k_{2}\hat{A}_{2}-6k_{2}\hat{A}_{3}-2k_{3}\hat{A}_{2}+\hat{A}_{1}\hat{d}_{0}-12k_{1}\hat{A}_{3},$
(38l) $\frac{\partial u_{2}}{\partial
C}=\frac{-1}{24(k_{1})^{2}}\frac{\partial k_{1}}{\partial
C}~{}~{},~{}~{}\frac{\partial u_{2}}{\partial
E}=\frac{-1}{24(k_{1})^{2}}\frac{\partial k_{1}}{\partial
E}~{}~{},~{}~{}\frac{\partial u_{2}}{\partial\kappa}=0,$ (38m) $\frac{\partial
v_{2}}{\partial C}=-2\left(k_{2}\frac{\partial\hat{A}_{2}}{\partial
C}+\hat{A}_{2}\frac{\partial k_{2}}{\partial
C}\right)-6\left(k_{2}\frac{\partial\hat{A}_{3}}{\partial
C}+\hat{A}_{3}\frac{\partial k_{2}}{\partial
C}\right)-2\left(k_{3}\frac{\partial\hat{A}_{2}}{\partial
C}+\hat{A}_{2}\frac{\partial k_{3}}{\partial C}\right)+\\\
\hat{A}_{1}\frac{\partial\hat{d}_{0}}{\partial
C}+\hat{d}_{0}\frac{\partial\hat{A}_{1}}{\partial
C}-12\left(k_{1}\frac{\partial\hat{A}_{3}}{\partial
C}+\hat{A}_{3}\frac{\partial k_{1}}{\partial C}\right),$ (38n) $\frac{\partial
v_{2}}{\partial E}=-2\left(k_{2}\frac{\partial\hat{A}_{2}}{\partial
E}+\hat{A}_{2}\frac{\partial k_{2}}{\partial
E}\right)-6\left(k_{2}\frac{\partial\hat{A}_{3}}{\partial
E}+\hat{A}_{3}\frac{\partial k_{2}}{\partial
E}\right)-2\left(k_{3}\frac{\partial\hat{A}_{2}}{\partial
E}+\hat{A}_{2}\frac{\partial k_{3}}{\partial E}\right)+\\\
\hat{A}_{1}\frac{\partial\hat{d}_{0}}{\partial
E}+\hat{d}_{0}\frac{\partial\hat{A}_{1}}{\partial
E}-12\left(k_{1}\frac{\partial\hat{A}_{3}}{\partial
E}+\hat{A}_{3}\frac{\partial k_{1}}{\partial E}\right),$ (38o) $\frac{\partial
v_{2}}{\partial\kappa}=-2k_{2}\frac{\partial\hat{A}_{2}}{\partial\kappa}-6k_{2}\frac{\partial\hat{A}_{3}}{\partial\kappa}-2k_{3}\frac{\partial\hat{A}_{2}}{\partial\kappa}+\hat{A}_{1}\frac{\partial\hat{d}_{0}}{\partial\kappa}-12k_{1}\frac{\partial\hat{A}_{3}}{\partial\kappa},$
(38p) $\frac{\partial\hat{d}_{0}}{\partial
C}=(2-\alpha)(\hat{c}_{0})^{1-\alpha}\frac{\partial\hat{c}_{0}}{\partial C},$
(38q) $\frac{\partial\hat{d}_{0}}{\partial
E}=(2-\alpha)(\hat{c}_{0})^{1-\alpha}\frac{\partial\hat{c}_{0}}{\partial
E}~{}~{},~{}~{}\frac{\partial\hat{d}_{0}}{\partial\kappa}=(2-\alpha)(\hat{c}_{0})^{1-\alpha}\frac{\partial\hat{c}_{0}}{\partial\kappa},$
(38r) $\frac{\partial\hat{c}_{0}}{\partial
C}=2\left(k_{4}\frac{\partial\hat{A}_{2}}{\partial
C}+\hat{A}_{2}\frac{\partial k_{4}}{\partial
C}\right)+k_{5}\frac{\partial\hat{A}_{1}}{\partial
C}+\hat{A}_{1}\frac{\partial k_{5}}{\partial C},$ (38s)
$\frac{\partial\hat{c}_{0}}{\partial
E}=2\left(k_{4}\frac{\partial\hat{A}_{2}}{\partial
E}+\hat{A}_{2}\frac{\partial k_{4}}{\partial
E}\right)+k_{5}\frac{\partial\hat{A}_{1}}{\partial
E}+\hat{A}_{1}\frac{\partial k_{5}}{\partial
E}~{}~{},~{}~{}\frac{\partial\hat{c}_{0}}{\partial\kappa}=2k_{4}\frac{\partial\hat{A}_{2}}{\partial\kappa}.$
(38t)

Next, we take derivative of (18jg) as follows

$\frac{\partial\hat{A}_{n+3}}{\partial C}=\frac{v_{3}\frac{\partial
u_{3}}{\partial C}-u_{3}\frac{\partial v_{3}}{\partial
C}}{(v_{3})^{2}}~{}~{},~{}~{}\frac{\partial\hat{A}_{n+3}}{\partial
E}=\frac{v_{3}\frac{\partial u_{3}}{\partial E}-u_{3}\frac{\partial
v_{3}}{\partial
E}}{(v_{3})^{2}}~{}~{},~{}~{}\frac{\partial\hat{A}_{n+3}}{\partial\kappa}=\frac{\frac{\partial
u_{3}}{\partial\kappa}}{v_{3}},$ (39a) where
$u_{3}=\left\\{-k_{2}(n+1)n-k_{3}(n+1)-k_{1}(n+1)n(n-1)\right\\}\hat{A}_{n+1}+\\\
\left\\{-k_{2}(n+2)(n+1)-2k_{1}(n+2)(n+1)n\right\\}\hat{A}_{n+2}+\hat{e}_{n},$
(39b) $v_{3}=k_{1}(n+3)(n+2)(n+1),$ (39c) $\frac{\partial u_{3}}{\partial
C}=\left\\{-k_{2}(n+1)n-k_{3}(n+1)-k_{1}(n+1)n(n-1)\right\\}\frac{\partial\hat{A}_{n+1}}{\partial
C}+\\\ \hat{A}_{n+1}\left\\{-\frac{\partial k_{2}}{\partial
C}(n+1)n-\frac{\partial k_{3}}{\partial C}(n+1)-\frac{\partial k_{1}}{\partial
C}(n+1)n(n-1)\right\\}+\\\
\left\\{-k_{2}(n+2)(n+1)-2k_{1}(n+2)(n+1)n\right\\}\frac{\partial\hat{A}_{n+2}}{\partial
C}+\\\ \hat{A}_{n+2}\left\\{-\frac{\partial k_{2}}{\partial
C}(n+2)(n+1)-2\frac{\partial k_{1}}{\partial
C}(n+2)(n+1)n\right\\}+\frac{\partial\hat{e}_{n}}{\partial C},$ (39d)
$\frac{\partial u_{3}}{\partial
E}=\left\\{-k_{2}(n+1)n-k_{3}(n+1)-k_{1}(n+1)n(n-1)\right\\}\frac{\partial\hat{A}_{n+1}}{\partial
E}+\\\ \hat{A}_{n+1}\left\\{-\frac{\partial k_{2}}{\partial
E}(n+1)n-\frac{\partial k_{3}}{\partial E}(n+1)-\frac{\partial k_{1}}{\partial
E}(n+1)n(n-1)\right\\}+\\\
\left\\{-k_{2}(n+2)(n+1)-2k_{1}(n+2)(n+1)n\right\\}\frac{\partial\hat{A}_{n+2}}{\partial
E}+\\\ \hat{A}_{n+2}\left\\{-\frac{\partial k_{2}}{\partial
E}(n+2)(n+1)-2\frac{\partial k_{1}}{\partial
E}(n+2)(n+1)n\right\\}+\frac{\partial\hat{e}_{n}}{\partial E},$ (39e)
$\frac{\partial
u_{3}}{\partial\kappa}=\left\\{-k_{2}(n+1)n-k_{3}(n+1)-k_{1}(n+1)n(n-1)\right\\}\frac{\partial\hat{A}_{n+1}}{\partial\kappa}+\\\
\left\\{-k_{2}(n+2)(n+1)-2k_{1}(n+2)(n+1)n\right\\}\frac{\partial\hat{A}_{n+2}}{\partial\kappa}+\frac{\partial\hat{e}_{n}}{\partial\kappa},$
(39f) $\frac{\partial v_{3}}{\partial C}=(n+3)(n+2)(n+1)\frac{\partial
k_{1}}{\partial C}~{}~{},~{}~{}\frac{\partial v_{3}}{\partial
E}=(n+3)(n+2)(n+1)\frac{\partial k_{1}}{\partial E}~{}~{},~{}~{}\frac{\partial
v_{3}}{\partial\kappa}=0,$ (39g) $\frac{\partial\hat{e}_{n}}{\partial
C}=\sum_{j=0}^{n}\hat{A}_{j}\frac{\partial\hat{d}_{n-j}}{\partial
C}+\hat{d}_{n-j}\frac{\partial\hat{A}_{j}}{\partial C},$ (39h)
$\frac{\partial\hat{e}_{n}}{\partial
E}=\sum_{j=0}^{n}\hat{A}_{j}\frac{\partial\hat{d}_{n-j}}{\partial
E}+\hat{d}_{n-j}\frac{\partial\hat{A}_{j}}{\partial E},$ (39i)
$\frac{\partial\hat{e}_{n}}{\partial\kappa}=\sum_{j=0}^{n}\hat{A}_{j}\frac{\partial\hat{d}_{n-j}}{\partial\kappa}+\hat{d}_{n-j}\frac{\partial\hat{A}_{j}}{\partial\kappa},$
(39j) $\frac{\partial\hat{d}_{n>0}}{\partial
C}=\frac{-1}{n(\hat{c}_{0})^{2}}\frac{\partial\hat{c}_{0}}{\partial
C}\sum_{j=1}^{n}(3j-\alpha j-n)\hat{c}_{j}\hat{d}_{n-j}+\\\
\frac{1}{n\hat{c}_{0}}\sum_{j=1}^{n}(3j-\alpha
j-n)\frac{\partial\hat{c}_{j}}{\partial
C}\hat{d}_{n-j}+\frac{1}{n\hat{c}_{0}}\sum_{j=1}^{n}(3j-\alpha
j-n)\hat{c}_{j}\frac{\partial\hat{d}_{n-j}}{\partial C},$ (39k)
$\frac{\partial\hat{d}_{n>0}}{\partial
E}=\frac{-1}{n(\hat{c}_{0})^{2}}\frac{\partial\hat{c}_{0}}{\partial
E}\sum_{j=1}^{n}(3j-\alpha j-n)\hat{c}_{j}\hat{d}_{n-j}+\\\
\frac{1}{n\hat{c}_{0}}\sum_{j=1}^{n}(3j-\alpha
j-n)\frac{\partial\hat{c}_{j}}{\partial
E}\hat{d}_{n-j}+\frac{1}{n\hat{c}_{0}}\sum_{j=1}^{n}(3j-\alpha
j-n)\hat{c}_{j}\frac{\partial\hat{d}_{n-j}}{\partial E},$ (39l)
$\frac{\partial\hat{d}_{n>0}}{\partial\kappa}=\frac{-1}{n(\hat{c}_{0})^{2}}\frac{\partial\hat{c}_{0}}{\partial\kappa}\sum_{j=1}^{n}(3j-\alpha
j-n)\hat{c}_{j}\hat{d}_{n-j}+\\\
\frac{1}{n\hat{c}_{0}}\sum_{j=1}^{n}(3j-\alpha
j-n)\frac{\partial\hat{c}_{j}}{\partial\kappa}\hat{d}_{n-j}+\frac{1}{n\hat{c}_{0}}\sum_{j=1}^{n}(3j-\alpha
j-n)\hat{c}_{j}\frac{\partial\hat{d}_{n-j}}{\partial\kappa},$ (39m)
$\frac{\partial\hat{c}_{n>0}}{\partial
C}=(k_{4}n+k_{5})(n+1)\frac{\partial\hat{A}_{n+1}}{\partial
C}+\left(\frac{\partial k_{4}}{\partial C}n+\frac{\partial k_{5}}{\partial
C}\right)(n+1)\hat{A}_{n+1}+\\\
k_{4}(n+2)(n+1)\frac{\partial\hat{A}_{n+2}}{\partial C}+\frac{\partial
k_{4}}{\partial C}(n+2)(n+1)\hat{A}_{n+2},$ (39n)
$\frac{\partial\hat{c}_{n>0}}{\partial
E}=(k_{4}n+k_{5})(n+1)\frac{\partial\hat{A}_{n+1}}{\partial
E}+\left(\frac{\partial k_{4}}{\partial E}n+\frac{\partial k_{5}}{\partial
E}\right)(n+1)\hat{A}_{n+1}+\\\
k_{4}(n+2)(n+1)\frac{\partial\hat{A}_{n+2}}{\partial E}+\frac{\partial
k_{4}}{\partial E}(n+2)(n+1)\hat{A}_{n+2},$ (39o)
$\frac{\partial\hat{c}_{n>0}}{\partial\kappa}=(k_{4}n+k_{5})(n+1)\frac{\partial\hat{A}_{n+1}}{\partial\kappa}+\left(\frac{\partial
k_{4}}{\partial\kappa}n+\frac{\partial
k_{5}}{\partial\kappa}\right)(n+1)\hat{A}_{n+1}+\\\
k_{4}(n+2)(n+1)\frac{\partial\hat{A}_{n+2}}{\partial\kappa}+\frac{\partial
k_{4}}{\partial\kappa}(n+2)(n+1)\hat{A}_{n+2}.$ (39p)

Finally, we take derivative of (16e) as follows

$\frac{\partial A_{1}}{\partial
C}=\frac{1}{1-\alpha}\left(\frac{u_{4}}{v_{4}}\right)^{\frac{\alpha}{1-\alpha}}\left(\frac{v_{4}\frac{\partial
u_{4}}{\partial C}-u_{4}\frac{\partial v_{4}}{\partial
C}}{(v_{4})^{2}}\right),$ (40a) $\frac{\partial A_{1}}{\partial
E}=\frac{1}{1-\alpha}\left(\frac{u_{4}}{v_{4}}\right)^{\frac{\alpha}{1-\alpha}}\left(\frac{v_{4}\frac{\partial
u_{4}}{\partial E}-u_{4}\frac{\partial v_{4}}{\partial
E}}{(v_{4})^{2}}\right)~{}~{},~{}~{}\frac{\partial A_{1}}{\partial\kappa}=0,$
(40b) where $u_{5}=k_{3}~{}~{},~{}~{}v_{5}=A_{0}(k_{5})^{2-\alpha},$ (40c)
$\frac{\partial u_{5}}{\partial C}=\frac{\partial k_{3}}{\partial
C}~{}~{},~{}~{}\frac{\partial u_{5}}{\partial E}=\frac{\partial
k_{3}}{\partial E}~{}~{},~{}~{}\frac{\partial u_{5}}{\partial\kappa}=0,$ (40d)
$\frac{\partial v_{5}}{\partial
C}=A_{0}(2-\alpha)(k_{5})^{1-\alpha}\frac{\partial k_{5}}{\partial
C}+(k_{5})^{2-\alpha}~{}~{},~{}~{}\frac{\partial v_{5}}{\partial
E}=A_{0}(2-\alpha)(k_{5})^{1-\alpha}\frac{\partial k_{5}}{\partial
E}~{}~{},~{}~{}\frac{\partial v_{5}}{\partial\kappa}=0.$ (40e)

## References

  * Sakiadis [1961] B. C. Sakiadis, “Boundary-layer behavior on continuous solid surfaces: II the boundary layer on a continuous flat surface,” AlChE J. 7, 221–225 (1961).
  * Weinstein and Ruschak [2004] S. J. Weinstein and K. J. Ruschak, “Coating flows,” Ann. Rev. Fluid Mech. 36, 29–53 (2004).
  * Blake, Clarke, and Ruschak [1994] T. D. Blake, A. Clarke, and K. J. Ruschak, “Hydrodynamic assist of dynamic wetting,” AIChE J. 40, 229–242 (1994).
  * Ganesh  _et al._ [2019] N. V. Ganesh, Q. M. Al-Mdallal, K. Reena, and S. Aman, “Blasius and sakiadis slip flow of $h_{2}o-c_{2}h_{6}o_{2}$ (50:50) based nanoliquid with different geometry of boehmite alumina nanoparticles,” Case Studies in Thermal Engineering 16 (2019), https://doi.org/10.1016/j.csite.2019.100546.
  * Azhar and Tasawar [2023] I. Azhar and A. Tasawar, “A study on heat transfer enhancement of copper ($cu$)-ethylene glycol based nanoparticle on radial stretching sheet,” Alexandria Engineering Journal 71, 13–20 (2023).
  * Abbas  _et al._ [2021] A. Abbas, I. Ijaz, M. Ashraf, and H. Ahmad, “Combined effects of variable density and thermal radiation on mhd sakiadis flow,” Case Studies in Thermal Engineering 28 (2021), https://doi.org/10.1016/j.csite.2021.101640.
  * Usafzai [2023] W. K. Usafzai, “Multiple exact solutions of second degree nanofluid slip flow and heat transport in porous medium,” Thermal Science and Engineering Progress 40 (2023), https://doi.org/10.1016/j.tsep.2023.101759.
  * Mishra  _et al._ [2023] P. Mishra, D. Kumar, Y. D. Reddy, and B. S. Goud, “Mhd williamson micropolar fluid flow pasting a non-linearly stretching sheet under the presence of non linear heat generation/ absorption,” Journal of the Indian Chemical Society 100 (2023), https://doi.org/10.1016/j.jics.2022.100845.
  * Usafzai and Aly [2023] W. K. Usafzai and E. H. Aly, “Multiple exact solutions for micropolar slip flow and heat transfer of a bidirectional moving plate,” Thermal Science and Engineering Progress 37 (2023), https://doi.org/10.1016/j.tsep.2022.101584.
  * Khazayinejad and Nourazar [2022] M. Khazayinejad and S. S. Nourazar, “On the effect of spatial fractional heat conduction in mhd boundary layer flow using $gr-fe_{3}o_{4}-h_{2}o$ hybrid nanofluid,” International Journal of Thermal Sciences 172 (2022), https://doi.org/10.1016/j.ijthermalsci.2021.107265.
  * Fox, Erickson, and Fan [1969] V. G. Fox, T. L. E. Erickson, and L. T. Fan, “The laminar boundary layer on a moving continuous flat sheet immersed in a non-newtonian fluid,” AlChE Journal 15, 327–333 (1969).
  * Blasius [1908] H. Blasius, “Grenzschichten in flussigkeiten mit kleiner reibung,” Zeitschrift fur Mathematik und Physik 56, 1–37 (1908).
  * Barlow  _et al._ [2017] N. S. Barlow, C. R. Stanton, N. Hill, S. J. Weinstein, and A. G. Cio, “On the summation of divergent, truncated, and underspecified power series via asymptotic approximants,” Q. J. Mech. Appl. Math. 70, 21–48 (2017).
  * Naghshineh  _et al._ [2023] N. Naghshineh, W. C. Reinberger, N. S. Barlow, M. A. Samaha, and S. J. Weinstein, “On the use of asymptotically motivated gauge functions to obtain convergent series solutions to nonlinear ODEs,” IMA Journal of Applied Mathematics (2023), doi:10.1093/imamat/hxad006.
  * Note [1] A gauge function is the usual independent variable in a series expansion. For example, in $e^{x}=\sum@\slimits@(x^{n}/n!)$, $x$ is the guage function. If we write $\textrm{exp}(e^{x})=\sum@\slimits@(e^{nx}/n!)$, $e^{x}$ is the gauge function [34, 35].
  * Bird, Armstrong, and Hassager [1987] R. B. Bird, R. Armstrong, and O. Hassager,  _Dynamics of Polymeric Liquids_ (John Wiley and Sons, 1987).
  * Schweizer [2022] P. Schweizer,  _Premetered Coating Methods: Attractiveness and Limitations_ (Springer, 2022).
  * Glass and Prud’homme [1997] E. Glass and R. K. Prud’homme, “In liquid film coating,” (S. F. Kistler and P. M. Schweitzer, New York: Chapman and Hall, 1997) Chap. Coating Rheology: Component Influence on the rheological response and performance in water-borne coatings in roll applications.
  * Pantokratoras [2015] A. Pantokratoras, “Non-similar blasius and sakiadis flow of a non-newtonian carreau fluid,” Journal of the Taiwan Institute of Chemical Engineers 56, 1–5 (2015).
  * Cebeci and Keller [1971] T. Cebeci and H. B. Keller, “Shooting and parallel shooting methods for solving the Falkner-Skan boundary-layer equation,” J. Comp. Phys. 7, 289–300 (1971).
  * Belden  _et al._ [2020] E. R. Belden, Z. A. Dickman, S. J. Weinstein, A. D. Archibee, E. Burroughs, and N. S. Barlow, “Asymptotic approximant for the falkner-skan boundary-layer equation,” Q. J. Mech. Appl. Math 73, 36–50 (2020).
  * Henrici [1956] P. Henrici, “Automatic computations with power series,” JACM 3, 10–15 (1956).
  * Churchill [1948] R. V. Churchill, “Complex variables,” (McGraw-Hill, 1948) Chap. VI: Power series.
  * Bataller [2010] C. R. Bataller, “Numerical comparisons of Blasius and Sakiadis flows,” MATEMATIKA 26, 187–196 (2010).
  * Fazio [2015] R. Fazio, “The iterative transformation method for the Sakiadis problem,” Comp. Fluids 106, 196–200 (2015).
  * Note [2] Some of the coefficients $a_{n}$ in (8d) are zero; hence, we use root test instead of ratio test.
  * Bender and Orszag [1978] C. M. Bender and S. A. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory_ (McGraw-Hill, 1978).
  * Pop and Gorla [1990] I. Pop and R. S. R. G. Gorla, “Second-order boundary layer solution for a continuous moving surface in a non-newtonian fluid,” Int. J. Engng Sci. 4, 313–322 (1990).
  * Van Dyke [1964] M. Van Dyke,  _Perturbation Methods in Fluid Mechanics_ (Academic, 1964).
  * Baker and Graves-Morris [1996] G. A. Baker and P. Graves-Morris,  _Padé Approximants_ (Cambridge, 1996).
  * Note [3] It is possible that poles of a Padé approximant arise within the physical domain of a problem for particular degrees of denominator and numerator. If the exact solution is expected to be finite within the physical domain, these Padés are deemed defective. For this exact reason, we did not use $M=25$ in Fig. 8 as it led to a defective approximant for the case of $\alpha=0.6$.
  * Note [4] For $\alpha=1$, one can use the Newtonian result (7c) in place of (16h); coefficients for (7c) are provided by Naghshineh et al[14].
  * Isaacson and Keller [1966] E. Isaacson and H. B. Keller,  _Analysis of Numerical Methods_ (John Wiley and Sons, New York, 1966).
  * Van Dyke [1975] M. Van Dyke, “Perturbation methods in fluid mechanics,” (Parabolic, 1975) Chap. 3.2: Gauge Functions and Order Symbols.
  * Leal [1992] L. G. Leal, “Laminar flow and convective transport processes. scaling principles and asymptotic analysis,” (Butterworth-Heinemann, 1992) Chap. 6.B: Asymptotic Expansions - General Considerations.