**3. HAM Solution**

Linear operators are

$$\text{L}\_f[f] = f'''' - f', \text{ L}\_{\mathfrak{F}}[\mathfrak{g}] = \mathfrak{g}'' - \mathfrak{g}, \text{ L}\_{\mathfrak{G}}[\mathfrak{g}] = \mathfrak{G}'' - \mathfrak{G}, \text{ L}\_{\mathfrak{G}}[\mathfrak{\phi}] = \mathfrak{\phi}'' - \mathfrak{\phi} \tag{24}$$

Initial guesses are taken as

$$f\_0[\zeta] = A \begin{bmatrix} 1 - e^{-\zeta} \end{bmatrix}, \ g\_0[\zeta] = \theta\_0[\zeta] = \phi\_0[\zeta] = e^{-\zeta} \tag{25}$$

with

$$\begin{aligned} \, \_f L\_f[f] \left[ p\_1 + p\_2 \epsilon^{-\zeta} + p\_3 \epsilon^{\zeta} \right] &= 0, \, \_L L\_{\xi}[g] \left[ p\_4 \epsilon^{-\zeta} + p\_5 \epsilon^{\zeta} \right] = 0, \\\, \_L L\_{\emptyset}[\Phi] \left[ p\_6 \epsilon^{-\zeta} + p\_7 \epsilon^{\zeta} \right] &= 0, \, \_L L\_{\Phi}[\Phi] \left[ p\_8 \epsilon^{-\zeta} + p\_9 \epsilon^{\zeta} \right] = 0 \end{aligned} \tag{26}$$

where *P*1, *P*2, *P*<sup>3</sup> ... , *P*<sup>9</sup> are called constants.

*Convergence Investigation*

By using the auxiliary parameters *f* , *<sup>g</sup>*, <sup>θ</sup> and <sup>φ</sup>, we analyzed the convergence regions for *f*'(ζ), *g*(ζ) and θ(ζ) of the modeled system of equations. At the 25th deformation order, the -− curves are presented here (see Figure 2). Convergence regions for *<sup>f</sup>* "(0), *<sup>g</sup>*'(0), <sup>θ</sup>'(0) and <sup>ϕ</sup>'(0) are <sup>−</sup>0.06 <sup>≤</sup> *<sup>f</sup>* ≤ 0.01, <sup>−</sup>0.08 <sup>≤</sup> *<sup>g</sup>* <sup>≤</sup> 0.03, <sup>−</sup>0.11 <sup>≤</sup> <sup>θ</sup> <sup>≤</sup> 0.05 and <sup>−</sup>0.11 <sup>≤</sup> <sup>φ</sup> ≤ 0.05, respectively.

**Figure 2.** -−curves for *f* "(0), *g*'(0), θ'(0) and ϕ'(0).
