Chebyshev Pseudospectral Method Overview

• The Chebyshev Pseudospectral Method (CPM) is a family of spectral techniques that approximates functions using Chebyshev polynomials and precise nodal grids.
• It employs both coefficient-space expansions and nodal collocation at Chebyshev–Gauss and Gauss–Lobatto nodes to discretize differential operators with spectral accuracy.
• CPM constructs almost banded system matrices to efficiently solve problems, addressing challenges like conditioning, singular sources, and fractional derivatives.

Chebyshev Pseudospectral Method (CPM) denotes a family of spectral discretizations in which functions are approximated by Chebyshev polynomials and differential operators are enforced either at selected collocation nodes or in coefficient space. In the cited literature, CPM encompasses Chebyshev–Gauss and Chebyshev–Gauss–Lobatto nodal schemes, shifted-Chebyshev operational-matrix constructions on finite intervals, time-space collocation, and the Chebyshev-coefficient “ultraspherical” method of Sheehan Olver and Alex Townsend for variable-coefficient ODEs (Olver et al., 2012). Across these variants, the common structure is the replacement of continuous derivatives, products, and boundary conditions by algebraic operators with spectral accuracy for sufficiently smooth data, together with problem-specific devices for conditioning, adaptivity, and the treatment of singular sources, nonhomogeneous data, fractional operators, and optimal-control constraints.

1. Basis functions, grids, and interval mappings

The foundational object is the Chebyshev polynomial of the first kind,

$$T_n(x)=\cos\!\bigl(n\arccos x\bigr),\qquad x\in[-1,1],$$

with the three-term recurrence

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$

Several formulations also use the second-kind family $U_n$, especially when derivatives or endpoint conditions are built directly into the basis [(Olver et al., 2012); (Ahrens et al., 26 May 2025)].

Two nodal sets dominate the literature. The Chebyshev–Gauss–Lobatto nodes are

$$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$

while the Chebyshev–Gauss nodes are

$$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$

The cited papers also use the sign-reversed Lobatto ordering

$$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$

which is the same grid with opposite orientation (Oltean et al., 2018). On $[0,1]$, the affine map

$$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$

produces Chebyshev–Gauss–Lobatto nodes adapted to shifted polynomials and diffusion problems on the unit square (Luchesi, 15 Sep 2025). On a general interval $[a,b]$, one uses

$x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$

and for moving time windows $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$0 the shifted basis is

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$1

with collocation instants obtained by mapping canonical Chebyshev nodes into physical time (Yousefian et al., 12 May 2025).

Two approximation viewpoints recur. In coefficient-space methods, one writes

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$2

so the unknown is the Chebyshev coefficient vector. In nodal collocation, one instead interpolates by Lagrange polynomials at Chebyshev nodes; the time-space Burgers–Fisher scheme, for example, uses the one-dimensional Chebyshev Lagrange basis

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$3

at Chebyshev–Gauss–Lobatto points in both $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$4 and $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$5 (Singh et al., 2023). The shifted-Chebyshev fractional formulation likewise expands

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$6

on $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$7 (Hoz et al., 14 Nov 2025).

2. Differentiation matrices, multiplication operators, and assembled systems

In nodal collocation, differentiation is realized by dense differentiation matrices. On a Chebyshev–Lobatto grid with $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$8 and $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$9 for $U_n$0, the first-derivative matrix has entries

$U_n$1

with diagonal entries defined by the row-sum identity

$U_n$2

Second derivatives are formed either as $U_n$3 or from the explicit Lobatto formula (Oltean et al., 2018). In time-space collocation, these one-dimensional operators are lifted to tensor form through Kronecker products: $U_n$4 (Singh et al., 2023).

The coefficient-space ultraspherical method uses a different algebra. If

$U_n$5

then multiplication by $U_n$6 becomes a Toeplitz $U_n$7 Hankel operator $U_n$8 of bandwidth $U_n$9. Higher derivatives are represented by diagonal operators $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$0 in the ultraspherical basis $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$1, and conversion operators $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$2 map $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$3-coefficients to $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$4-coefficients. The $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$5th term of a variable-coefficient ODE therefore acts as

$$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$6

on the Chebyshev-coefficient vector (Olver et al., 2012).

A central structural consequence is “almost bandedness.” After imposing $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$7 boundary conditions by row replacement and permuting them to the top, the global matrix contains $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$8 dense boundary rows, a fixed number of diagonals from differentiation and conversion, and a multiplication band of width $$x_i=\cos\!\Bigl(\frac{i\pi}{N}\Bigr),\qquad i=0,\dots,N,$$9. Because $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$0 and the differential order $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$1 are fixed by the smoothness of the coefficients, the discretized system is almost banded rather than fully dense (Olver et al., 2012).

The same operator logic extends to higher dimensions and mixed bases. In the transient anisotropic diffusion formulation on $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$2, the full Chebyshev differentiation matrix $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$3 is sliced into interior-only blocks $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$4 and $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$5, and anisotropy enters through diagonal matrices of sampled conductivities together with boundary corrections derived from Robin-type conditions (Luchesi, 15 Sep 2025). In superconducting strip problems on $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$6, $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$7, the method combines Chebyshev differentiation in the finite $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$8-direction with Hermite-function differentiation in the unbounded $$x_i=\cos\!\Bigl(\frac{(2i+1)\pi}{2N+2}\Bigr),\qquad i=0,\dots,N.$$9-direction, producing

$$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$0

for the two coordinate directions (Sokolovsky et al., 2021).

3. Conditioning, preconditioning, and solver architecture

The best-developed conditioning theory in the cited set is the Olver–Townsend ultraspherical scheme for linear ODEs with variable coefficients. The key diagonal preconditioner is

$$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$1

which balances the growth of ultraspherical differentiation and conversion so that

$$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$2

with $$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$3 compact in an appropriate $$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$4 norm (Olver et al., 2012). For Dirichlet boundary conditions, one may take $$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$5, yielding a bounded $$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$6-norm condition number of $$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$7 uniformly in $$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$8. In that sense, the method is well conditioned in the standard Euclidean norm for that class of boundary data.

The associated direct solver is an adaptive QR factorization applied to an “infinite” almost-banded operator. Givens rotations are introduced column by column, each filled-in row is maintained in $$x^{(s)}_j=-\cos\!\Bigl(\frac{\pi j}{N}\Bigr),\qquad j=0,\dots,N,$$9 storage, and truncation is determined by the transformed right-hand side. The stopping test is

$[0,1]$0

with a representative tolerance $[0,1]$1. The operation count is $[0,1]$2, the memory requirement is $[0,1]$3, and the reported implementation can “efficiently and reliably solve for solutions that require as many as a million unknowns” (Olver et al., 2012).

Other CPM variants emphasize stability through different mechanisms. The time-space Burgers–Fisher scheme derives a discrete weighted-norm estimate,

$[0,1]$4

after discrete integration by parts, thereby proving spectral stability of the fully coupled time-space collocation method (Singh et al., 2023). The shifted-Chebyshev fractional construction identifies the dominant numerical hazard not in the discrete-cosine transform, whose conditioning is $[0,1]$5, but in the coefficient matrix $[0,1]$6 used to generate operational matrices. Because its entries grow like $[0,1]$7 and alternate in sign, the method uses symbolic or variable-precision arithmetic with digits chosen approximately linearly in $[0,1]$8 to preserve stability up to $[0,1]$9 in the reported tests (Hoz et al., 14 Nov 2025).

4. Boundary conditions, singular sources, and PDE generalizations

A persistent theme in CPM research is that singular or nonhomogeneous data are typically not inserted naively into a single global discretization. The “Particle-without-Particle” method is a direct illustration. For a one-dimensional linear PDE with a source

$$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$0

the domain is split into $$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$1 and $$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$2, each with its own Chebyshev–Lobatto grid. One solves the homogeneous PDE $$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$3 in each subdomain and replaces the singular source by jump conditions at the interface. For the model problem

$$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$4

the jumps are

$$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$5

(Oltean et al., 2018). The paper proves that the method applies to any linear PDE whose source is any linear combination of delta distributions and derivatives thereof supported on a one-dimensional subspace of the problem domain.

For nonlinear initial-boundary-value problems, nonhomogeneous data are often homogenized before collocation. In the time-space Chebyshev pseudospectral method for the generalized Burgers–Fisher equation,

$$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$6

a lifting function $$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$7 interpolates the initial and boundary data, and the unknown is rewritten as

$$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$8

so that $$x_i=\frac12\Bigl[1-\cos\!\Bigl(\frac{i\pi}{n}\Bigr)\Bigr]$$9 satisfies homogeneous initial and boundary conditions (Singh et al., 2023). Interior collocation then yields a nonlinear algebraic system of dimension $[a,b]$0, solved by Newton–Raphson.

In transient anisotropic diffusion, the spatial operator is built from a time-dependent rotated diffusion tensor

$[a,b]$1

with $[a,b]$2 determined by a rotation angle $[a,b]$3. Robin-type conditions on all four sides of the square are reformulated to express normal derivatives in terms of $[a,b]$4 and boundary data, after which lexicographic ordering yields the semidiscrete system

$[a,b]$5

Time advancement is performed by Crank–Nicolson, requiring a solve with $[a,b]$6 at each step when $[a,b]$7 varies in time (Luchesi, 15 Sep 2025).

The Hermite–Chebyshev superconducting-strip formulation shows a further geometric generalization. On $[a,b]$8, $[a,b]$9, the unknown is expanded in a tensor basis of Hermite functions and Chebyshev polynomials, while a stream function $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$0 is expanded in second-kind Chebyshev polynomials $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$1 so that the Dirichlet condition $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$2 is built into the basis. At each time step, the method computes current density from spectral derivatives, evaluates a nonlinear current-voltage law pointwise, solves $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$3 dense systems of size $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$4, reconstructs the field, and advances with an ODE solver such as MATLAB’s ode15s (Sokolovsky et al., 2021).

5. Fractional operators, optimization, control, and online identification

The CPM framework has also been extended well beyond classical integer-order boundary-value problems. For Caputo-type advection-diffusion equations, shifted Chebyshev polynomials on $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$5 are used to construct an operational matrix $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$6 for the Caputo derivative and an operational matrix $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$7 for the Riemann–Liouville integral. If

$x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$8

then at the shifted Chebyshev–Gauss–Lobatto nodes $x_j=\frac12(a+b)+\frac12(b-a)\,x^{(s)}_j,$9,

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$00

For the one-dimensional PDE

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$01

the discretization leads to a Sylvester equation

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$02

solved by MATLAB’s built-in lyap; in $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$03 spatial dimensions it becomes a tensor Sylvester equation solved by sylvesterND from Cuesta–de la Hoz ’24 (Hoz et al., 14 Nov 2025).

In direct optimal control, the cited methodology transcribes an ODE-constrained problem into a nonlinear program by approximating only the highest-order derivative in a Chebyshev series and recovering lower-order states by successive integration. If $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$04 is approximated as

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$05

then sparse integration matrices $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$06 map the coefficient vector $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$07 to lower derivatives and states. Dynamic constraints are collocated at Chebyshev nodes, endpoint constraints are enforced directly for CGL or by interpolation for CG, and the cost is integrated by Clenshaw–Curtis weights. The decision vector may take the form

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$08

(Ahrens et al., 26 May 2025).

A related but distinct line of work uses CPM for one-dimensional optimization. The exact line-search method of “Optimization via Chebyshev Polynomials” constructs a fourth-order Chebyshev interpolant on five CGL points, differentiates it to obtain a cubic model, classifies its roots analytically, and then refines the candidate minimizer by a Newton step using first- and second-order Chebyshev pseudospectral differentiation matrices. If the Hessian surrogate is non-positive or the iterate exits the interval, the method falls back to a golden-section or Brent step. The paper states that the resulting update converges quadratically when second-order information is used (Elgindy, 2016).

The online identification work extends the method to moving time windows. Over $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$09, the system drift is approximated as

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$10

with $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$11 Chebyshev sampling instants. Coefficients are computed by the regularized least-squares formula

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$12

and continuity across windows is imposed by matching all derivatives up to order $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$13 at $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$14. The same piecewise-Chebyshev approximation is then embedded in an observer whose gain solves a Lyapunov equation, and the paper proves uniform boundedness of the parameter and state errors (Yousefian et al., 12 May 2025).

6. Convergence behavior, quantitative results, and methodological scope

The dominant accuracy claim across the cited literature is spectral or exponential convergence for sufficiently smooth or analytic solutions. In the Particle-without-Particle formulation, if each subdomain solution $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$15 is smooth, then Chebyshev collocation yields

$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$16

for some $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$17 determined by analyticity-strip width (Oltean et al., 2018). The time-space Burgers–Fisher study reports max-norm errors decaying from $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$18 at $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$19 to $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$20 at $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$21, which is the expected spectral decay until round-off (Singh et al., 2023). The anisotropic diffusion abstract states that analytic solutions exhibit $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$22 accuracy, whereas solutions with only $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$23 continuous derivatives recover $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$24 algebraic decay (Luchesi, 15 Sep 2025). The Hermite–Chebyshev strip method is spectrally convergent in $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$25 but only algebraically convergent in $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$26 when the solution decays with a power law at infinity (Sokolovsky et al., 2021).

Several papers report concrete performance data.

Context	Reported outcome	Source
Variable-coefficient ODE $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$27, $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$28	adapts to $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$29 for 16-digit accuracy; condition number of $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$30 stays $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$31; solve time $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$32	(Olver et al., 2012)
Airy BVP with $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$33	solver finds $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$34 in $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$35 (C++), with uniform error $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$36	(Olver et al., 2012)
Advection equation with moving $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$37 source	at $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$38, PwP + Chebyshev PSC gives $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$39 error $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$40 versus $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$41 for $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$42-approximation + 6th-order finite difference	(Oltean et al., 2018)
Elliptic Poisson problem with ring $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$43 source	at $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$44, PwP + CL-PSC gives $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$45 error $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$46 versus $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$47 for Tornberg–Kreiss $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$48-approximation	(Oltean et al., 2018)
Inhomogeneous superconducting strip	$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$49 gives $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$50; $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$51 reduces error to $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$52	(Sokolovsky et al., 2021)
Superconducting dynamo benchmark	voltage error $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$53 with $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$54, CPU time $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$55–$$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$56 min/cycle	(Sokolovsky et al., 2021)
Rocket landing flip maneuver	with $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$57, Chebyshev CG/CGL/second-kind and Legendre LGL give final mass $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$58 kg in $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$59 s; differences $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$60 kg	(Ahrens et al., 26 May 2025)
One-dimensional optimization tests	CPSLSM achieved 10–15 correct digits versus MATLAB’s fminbnd typically 5–8 digits; iterations 2–15 versus 5–40	(Elgindy, 2016)

These results also clarify the scope of the method. The cited literature documents several counterexamples to the view that CPM is limited to a single dense differentiation matrix on one interval. The method appears as an almost-banded coefficient-space scheme for linear ODEs (Olver et al., 2012), a multi-domain collocation method with exact jump conditions for distributional sources (Oltean et al., 2018), a fully coupled time-space collocation method (Singh et al., 2023), a shifted-Chebyshev operational-matrix construction for fractional derivatives (Hoz et al., 14 Nov 2025), a mixed Hermite–Chebyshev method on $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$61 (Sokolovsky et al., 2021), and an integral-collocation transcription for nonlinear programming in optimal control (Ahrens et al., 26 May 2025).

At the same time, the literature records clear limitations. If coefficient functions develop singularities near $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$62, the number $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$63 of modes required to represent them may grow with $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$64, and the direct ultraspherical solver remains stable but its cost becomes $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$65 (Olver et al., 2012). Highly oscillatory temporal data in the fractional setting require $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$66 to resolve frequencies $$T_0(x)=1,\qquad T_1(x)=x,\qquad T_{n+1}(x)=2x\,T_n(x)-T_{n-1}(x).$$67 (Hoz et al., 14 Nov 2025). In infinite-domain strip problems, the Hermite direction loses spectral convergence when the underlying solution has only algebraic decay (Sokolovsky et al., 2021). These caveats do not invalidate CPM; rather, they delimit the regimes in which its characteristic spectral efficiency is fully realized.