In this section, we describe the algorithm to describe how the multidomain bivariate spectral collocation method can be applied to solve nonlinear parabolic PDEs. We shall consider a general second-order nonlinear PDE,

### 2.1. The quasi-linearization method

The quasi-linearization method (QLM) of Bellman and Kalaba [18] is a technique that is used to simplify nonlinear ordinary and partial differential equations. The technique has been adopted and generalized in further studies presented in [22, 23]. The QLM is based on the Newton-Raphson method and is constructed from the linear terms of Taylor series expansion about an initial approximation to solution. The QLM assumes that the difference between solutions at two successive iterations denoted by *u*_{s} and *u*_{s + 1} is very small. Applying the QLM on Eq. (1) yields

where prime denotes differentiation with respect to *x* and *s* denotes the iteration level. Eq. (4) can be written in compact form as

where *u*^{(0)} = *u*. Using the expanded form of Eq. (5) in Eq. (1), we obtain the QLM scheme for approximating the solution *u*_{s + 1}(*x*, *t*) at the (*s* + 1)th iteration level as

where

The dot here denotes differentiation with respect to the time *t*. Starting with an initial approximation *u*_{0}, the QLM scheme is solved iteratively until a solution with desired accuracy requirements is obtained. The multidomain approach is implemented on the linearized scheme (6) as illustrated below. For the purpose of this study, we shall apply the multidomain approach on the time (*t*) variable only.

Let *t* ∈ *Γ* where *Γ* ∈ [0, *T*]. The domain *Γ* is decomposed into *p* nonoverlapping intervals as

The domain *t* ∈ [*t*_{k − 1}, *t*_{k}] in each of the *k*th subdomain is first transformed to *τ* ∈ [−1, 1] using the linear transformation

before the spectral collocation is applied. Similarly, the spatial domain *x* ∈ [*a*, *b*] is transformed to *η* ∈ [−1, 1] using the linear transformation

The collocation nodes are the symmetrically distributed Gauss-Lobatto grid points defined on the interval [−1, 1] by,

To distinguish between the solutions at different subdomains we shall use, *u*^{(k)}, *k* = 1, 2, …, *p*, to denote solution at the *k*th subinterval. The PDE is solved independently in each subinterval. In the first subinterval we must solve,

subject to boundary and initial conditions

After the solution in the first interval *Γ*_{1} has been computed, the solutions at the subsequent *k*th subinterval are computed by using the solution at the right hand boundary of the (*k* − 1)*th* interval as an initial solution. Thus in the next subintervals, *k* = 2, 3, …, *p*, we must solve

subject to boundary and initial conditions

In the solution process, the approximate solution that is searched for takes a form of a bivariate Lagrange interpolation polynomial. The solution at each subinterval is approximated as

The first and second spatial derivatives are evaluated at the collocation nodes (*η*_{i}, *τ*_{j}) for *j* = 0, 1, 2, …, *M* as follows

where D^=b−a2Dof size (*N* + 1) × (*N* + 1) is the standard first-order Chebyshev differentiation matrix as defined in [15]. The time derivative is evaluated at the collocation nodes (*η*_{i}, *τ*_{j}) for *i* = 0, 1, 2, …, *N* as

where d^j,q=tk−tk−12dj,q, *j*, *q* = 0, 1, 2, …, *M* of size (*M* + 1) × (*M* + 1) is the standard first-order Chebyshev differentiation matrix,

and *T* denotes matrix transpose. Using the definitions (17)–(18), we express Eq. (6) in matrix form as

By changing the indices, Eq. (20) can be written as

where

Eq. (21) constitutes an *M*(*N* + 1) × *M*(*N* + 1) matrix system given by

where

and **I** is an identity matrix of size (*N* + 1) × (*N* + 1). The boundary conditions at the collocation points are

These boundary conditions are imposed on the main diagonal submatrices of the matrix system (23) to obtain a new system which takes the form

where

The matrix system (26) is solved for *U*^{(k)}, *k* = 1, 2, …, *p*. The solutions at different subdomains are matched together along common boundaries to give the desired approximate solution. The patching condition is given by

which denotes the solution at the boundaries of the subintervals.