Boundary-eliminated Chebyshev collocation for a two-point Dirichlet BVP

Input an integer $N\ge 2$, real functions $a,b,f$ given at the Chebyshev–Gauss–Lobatto nodes $x_j=\cos (j\pi \mathrm{/}N)$, and endpoint data $g_{-},g_{+}\in \mathbb{R}$. The continuous BVP being discretized is

$$-u^{\mathrm{\prime }\mathrm{\prime }}(x)+a(x)u^{\mathrm{\prime }}(x)+b(x)u(x)=f(x),-1<x<1,u(-1)=g_{-},u(1)=g_{+}\mathrm{.}$$

The finite algorithm uses only the values at the nodes and does not assume a continuous BVP solvability theorem. Construct the matrix used internally by Chebyshev–Gauss–Lobatto polynomial differentiation matrix explicitly: form $L_j(x)={\prod }_{0\le k\le N,k\mathrm{\ne }j}(x-x_k)\mathrm{/}(x_j-x_k)$ and set $(D_N{)}_{ij}=L_j^{\mathrm{\prime }}(x_i)$ for $0\le i,j\le N$. Then set $D_N^{(2)}=D_N{D}_N$. For each interior row $i=1,\dots ,N-1$ and each column $j=0,\dots ,N$, put

$$C_{ij}=-(D_N^{(2)}{)}_{ij}+a(x_i)(D_N{)}_{ij}+b(x_i){\delta }_{ij}\mathrm{.}$$

Let $A=(C_{ij}{)}_{1\le i,j\le N-1}$ and $r_i=f(x_i)-C_{i0}{g}_{+}-C_{iN}{g}_{-}$. If $A$ is singular, return status singular-interior-matrix and the matrix $A$. If $A$ is invertible, solve $A{U}_I=r$, set $U_0=g_{+}$ and $U_N=g_{-}$, and return success with $U=(U_0,U_1,\dots ,U_N)$. This exact-arithmetic finite solve asserts no unconditional nonsingularity, stability, convergence, or PDE error bound.