Mathematically, we will have to solve the following nonlinear boundary value problem:

where $𝒩$ is a third-order nonlinear differential operator of the form,

which may also possibly depend on some real parameters. The right-hand side member $f$ is a continuous function of all its arguments, and $A,$ $B$, and $C$ are specified constants.

Each smooth solution to the problem (1) is mapped to the interval $[-1,1]$ and represented by an expansion in Chebyshev polynomials of the form

where $N$ is a natural number, the so called the resolution, and $T_j\left(x\right):=cos(j\ast arccos(x))$ are the well known Chebyshev polynomials.

When constructing a solution (Chebfun object), the system computes the coefficients $u_j$ by interpolating the target function $u_{Chebfun}$ at $N+1$ Chebyshev’s extreme points,

The polynomial degree $N$ is automatically determined so that the above representation is as accurate as possible in double-precision arithmetic.

Polynomial interpolation in Chebyshev nodes is known to be near-optimal for approximating functions that are smooth, converging geometrically for analytic functions. The Fast Fourier transforms (FFTs) can be used to map function values $u_{Chebfun}(x_j)$ to coefficients $u_j$ and vice versa, in $O(NlogN)$ operations.

The construction process actually begins by sampling the target function at $2^n+1$ points, with $n=3,4,\mathrm{\dots }$. The optimal degree $N$ is then determined such that $|u_j|$ is close to zero, relative to the coefficient of largest magnitude, for all $j>N.$

Chebfun is quite thoroughly described in the monographs of L. N. Trefethen and his students [8] and [9] (see also http://www.chebfun.org ). The motto under which this programming environment is further built is ”feel symbolic but run at the speed of numerics”. The aim of this open-source object-oriented MATLAB package is to achieve for functions what floating-point arithmetic achieves for numbers.

The fundamental principle of the Chebfun system is to evaluate functions in sufficiently many Chebyshev nodes (3) for a polynomial interpolant to be accurate to machine precision. Whenever the function of interest is only defined implicitly, as the solution of a differential equation, Chebfun works similarly with MATLAB based on collocation in the Chebyshev nodes and lazy evaluation of the associated spectral differentiation matrices.

This programming environment has proven to be very performing but sometimes encounters difficulties when the analyzed non-linearity contains rational (fractional) powers. To overcome this inconvenience, we initiated a method of continuation, that is, we solved simpler samples that do not contain such powers, whose solutions then became initial data for more complex problems.

To prove the convergence of the solution it is very valuable the estimation

which gives us a measure of how accurately Chebfun performed. Final error estimates for differential equations and boundary conditions are provided separately. They are generally better than the error defined by (4).

The error introduced with the formula (4) will help us to find an optimal length of the integration interval. This will be done for each example below without explicitly mentioning it.

To have precise definitions for classifying the rate of convergence of Chebfun in solving various problems, we refer to the well-known monograph of [10] (Section 2.3). Roughly, an expansion (series) of type (2) is said to have the property of supergeometric, geometric, or subgeometric convergence depending upon whether;

In a log-linear plot, $log\left|u_n\right|$ vs.$n$, it is relatively simple to distinguish between these three types of convergence. We will do this for each of the examples below.

A very useful example of Chebfun’s code is found in our contribution [11]. There we solved a singular nonlinear BVP of the third order that accepts a prime integral and comes from the simulation of the large-scale circulation in an ocean (the Gulf Stream).

Let us consider the Blasius equation

along with the impermeability of the wall, the nonslip and the free stream velocity outside the boundary layer, respectively, boundary conditions

where the primes denote differentiation concerning the independent variable $x.$

Here we merely note that the function $u^{\prime }\left(x\right)$ represents, after suitable normalization, the fluid velocity parallel to the wall (the Ox axis) and $u$ stands for normalized stream function. This problem is frequently dealt with in fluid dynamics textbooks but very rarely solved as a BVP.

Roughly, this problem describes the steady two-dimensional laminar boundary layer that forms on a semi-infinite plate which is held parallel to a constant unidirectional flow.

As a member of the Chebfun team, the author of [12] shows that the solution to this problem is a smooth, monotonically increasing function that rapidly converges to a linear polynomial away from the origin. The singularity of this solution is established at $x=-5.6900380545$.

Our Chebfun solution along with its first two derivatives is depicted in Fig. 1. From the left panel of Fig.2 we infer that Newton’s method has second order and from the right panel of the same figure we can state that the convergence of spectral method tends to a geometric one. How the coefficients decrease is a correct indication of the convergence of the solution.

However, from the fifth line of Table 1 we see that the numerical value of the second derivative of solution in origin, the so-called the local shear stress, reaches the exact value computed in [13].

We now attach to the Blasius equation (6) the following slip boundary condition

where $c$ is a non-dimensional parameter.

In our computations, we have used $c:=2.0$ and obtained local shear stress $u^{\prime \prime }\left(0\right)=2.396618681578377e-01.$ Chebfun needed a modest resolution $N=33$ and provided the solution with an accuracy of order $8.568214448343653e-09$ in $8$ iterations.

The solution to this problem is depicted in Fig. 3 and the convergence is argued in Fig. 4. It is visible from the second figure that the Newton method is second-order and the convergence of the solution is clearly geometric.

The length of the integration domain has been $X:=6$ but for $X$ in the range $[4,20]$ we have obtained almost identical numerical results, i.e., that differs by less than ${10}^{-6}.$

It is really surprising to note that the behaviour of a certain type of electromechanical energy converter (EMEC) is governed by the same Blasius equation (6) (see [14]). However, the behaviour of the solutions to the EMEC problem is essentially different from the boundary layer solution because of different initial and boundary conditions. It is not our purpose to analyze IVPs so we will not dwell on this problem.

The laminar boundary-layer behavior on a moving continuous flat surface is now investigated. We are concerned now with the classical Blasius equation (6) for the flat plate of infinite length along with the boundary conditions

(see also the paper of Sakiadis [15]).

For this problem, Chebfun worked with $N=33$ in (2) in order to get an accuracy of order $2.711625848703521e-10$ and Newton’s method has carried out the solution in $7$ iterations. The local shearing stress $u^{\prime \prime }\left(0\right)$ obtained, equals the value $-4.439043651681859e-01.$

The Chebfun solution for this problem along with its first two derivatives is depicted in Fig. 5. From the left panel of Fig. 6 we observe that Newton’s method has an order less than two but close to this value and from the right panel of the same figure we can state that the convergence of spectral method tends to a geometric rate.

The well-known Falkner-Skan problem reads

where $\lambda$ is a real parameter (wedge angle). It generalizes Blasius’ flow to wedge flow, i.e. flows in which the plate is not parallel to the flow.

The existence of a solution to the problem (10) when $\lambda >0$ was first proved by Weyl in [16] and the uniqueness for $0\le \lambda \le 1$ is analyzed in Craven and Peletier, see [17].

For this problem, Chebfun worked with $N=26$ in (2) in order to get an accuracy of order $6.991519876438462e-09$ and Newton’s method has carried out the solution in $6$ iterations with a convergence rate of order $2.$

The assumed value for $\lambda$ has been $0.65$ and the local shearing stress obtained equals the value

There are two important sources by which we can verify this value. The first is due to [18] who shows that this value must be between $4\lambda /3$ and $4\left(\lambda +1\right)/3,$ which happens for the value found by Chebfun.

Moreover, by a simple integration of the Falkner-Skan equation and the introduction of the boundary conditions from (10), the following important relationship is found in [19]

Chebfun calculates the two integrals simply and elegantly, and the value found for shear stress with (12) coincides with the exact one (11) up to the sixth decimal place.

The solution to the problem (10) is depicted in Fig. 7 and details about the convergence of Newton’s method and the solution can be found in Fig. 8. It is clear that Newton’s method is of second order and the convergence of the solution approaches the geometric order.

In this section, we will be concerned with the problem classified as difficult, namely

with $n$ a sub-unit parameter. In his paper [20] Liao considers that this problem is a challenging nonlinear one for numerical techniques.

We want to note from the beginning that, a direct approach to this issue using Chebfun, led to a failure. Chebfun warns that the Newton method failed and asks for a better initial guess.

In this situation, we resorted to a kind of continuation, i.e. we first solved the problem with $n:=1,$ i.e. Blasius problem, and then we used its solution as an initial guess for the problem with an arbitrary $n$.

The results are presented in Table 1 and a set of five solutions are presented in Fig. 9. It is observed that the error in the calculation of the solutions is of an order of ${10}^{-8}$ or better. The order of the Newtonian method has not been reported, but it tends towards two in every case.

The numerical results have been stable for lengths $X\in [8,40].$

Along with Magyari and Keller [21], we consider the singular nonlinear BVP

where $u_w=0$ corresponds to an impermeable wall, $u_w>0$ to suction and to lateral injection of the fluid through a permeable wall. The coefficient $m$ comes from the power-law variations of the stretching velocity.

In our numerical experiments, the assumed value for $m,$ has been $-1/3.$ For this value, in [21] the authors find an exact solution with the maximum value of its first derivative of $1+u_w^2/2.$ Our numerical solutions reported in Fig. 10 visibly confirm this estimation.

Heidel in [22] observes that equation (14) differs only slightly from the Falkner-Skan equation (10). However, a basic difference between (14) and the Falkner-Skan equation is that the former is invariant under the transformation $y\left(x\right):=cu\left(cx\right)$ whereas the Falkner-Skan equation is not. This allows the reduction order of (14).

A fully analogous problem from [23] has been solved. To save space, we do not present these results.

King and Lewellen consider in [25] the following singular BVP, i.e. the nonlinear differential system

along with boundary conditions

where $n$ is the power of the radius in the asymptotic assumption of external axial velocity $\left(-1\le n\le 1\right)$ and $S$ is a physical parameter proportional with the intensity of the magnetic field ( $S=O\left(1\right)$ is assumed).

Physically, equations (16), with the boundary conditions of equation (17), determine the velocities in the boundary layer as a function of the similarity variable $x$ and the two parameters $n$ and $S$.

Our numerical results are depicted in Fig. 11. They have been obtained for $X:=20$ but are perfectly stable for $X\in [10,40].$ In fact, the value $X=:40$ was also used in the work [24]. Unfortunately, a direct and quantitative comparison with these results is almost impossible to track. However, the allure of the variation of the unknowns $u$ and $v$ are the same.

It is interesting to note that an elementary solution of the system is given by the pair $(0,1)$, and the parameter $S$ affects only the linear part of this system. Our numerical experiments confirm those from the paper [25] which observe that for $S\ne 0$, the oscillations of $u$ and $v$ are damped out and the thickness of the boundary layer drastically decreases.