*3.1. Procedure for Solving the Direct TVIDE Problem*

First, we linearize (8) by uniformly discretizing the temporal domain into *M* subintervals with time step *<sup>τ</sup>*. Then, we specify (8) at a time *tm* <sup>=</sup> *<sup>m</sup><sup>τ</sup>* for *<sup>m</sup>* <sup>∈</sup> <sup>N</sup> and use the first-order forward difference quotient to estimate the time derivative term *ut*. Next, we replace each *x* by *xk* for *k* ∈ {1, 2, 3, ..., *N*} as generated by the zeros of the shifted Chebyshev polynomial *SN*(*x*) defined in (2). Thus, we have

$$\frac{\mu^{\langle m\rangle} - \mu^{\langle m-1\rangle}}{\tau} + \mathcal{L}u^{\langle m\rangle} = \int\_0^{t\_m} \mathbf{x}\_1(\mathbf{x}\_k, \eta) u(\mathbf{x}\_k, \eta) d\eta + \int\_0^{\mathbf{x}\_k} \mathbf{x}\_2(\xi, t\_m) u(\xi, t\_m) d\xi + F^{\langle m\rangle},\tag{11}$$

where *um* = *um*(*xk*) = *u*(*xk*, *tm*) and *Fm* = *Fm*(*xk*) = *F*(*xk*, *tm*). Next, consider the first integral term with respect to time by letting it be *J m* <sup>1</sup> (*xk*), we approximate *J m* <sup>1</sup> (*xk*) by using the trapezoidal rule. Thus, we approximate *J m* <sup>1</sup> (*xk*) as

$$\begin{split} I\_{1}^{\{m\}}(\mathbf{x}\_{k}) &:= \int\_{0}^{t\_{m}} \kappa\_{1}(\mathbf{x}\_{k},\eta) u(\mathbf{x}\_{k},\eta) d\eta \\ &= \sum\_{i=0}^{m-1} \int\_{t\_{i}}^{t\_{i+1}} \kappa\_{1}(\mathbf{x}\_{k},\eta) u(\mathbf{x}\_{k},\eta) d\eta \\ &\approx \sum\_{i=0}^{m-1} \frac{\tau}{2} \left( \kappa\_{1}^{(i)}(\mathbf{x}\_{k}) u^{(i)}(\mathbf{x}\_{k}) + \kappa\_{1}^{(i+1)}(\mathbf{x}\_{k}) u^{(i+1)}(\mathbf{x}\_{k}) \right) \\ &= \frac{\tau}{2} \kappa\_{1}^{(0)}(\mathbf{x}\_{k}) u^{(0)}(\mathbf{x}\_{k}) + \tau \sum\_{i=1}^{m-1} \kappa\_{1}^{(i)}(\mathbf{x}\_{k}) u^{(i)}(\mathbf{x}\_{k}) + \frac{\tau}{2} \kappa\_{1}^{(m)}(\mathbf{x}\_{k}) u^{(m)}(\mathbf{x}\_{k}) \end{split}$$

for each *xk* ∈ {*x*1, *x*2, *x*3, ..., *xN*}. The above equation can be written, in matrix form, as

$$\mathbf{J}\_1^{\langle m\rangle} = \frac{\tau}{2} \mathbf{K}\_1^{\langle 0\rangle} \mathbf{u}^{\langle 0\rangle} + \tau \sum\_{i=1}^{m-1} \mathbf{K}\_1^{\langle i\rangle} \mathbf{u}^{\langle i\rangle} + \frac{\tau}{2} \mathbf{K}\_1^{\langle m\rangle} \mathbf{u}^{\langle m\rangle},\tag{12}$$

where each parameter in (12) can be defined as follows:

$$\begin{array}{lcl} \mathbf{J}\_{1}^{\langle m \rangle} &=& \left[ l\_{1}^{\langle m \rangle}(\mathbf{x}\_{1}), l\_{1}^{\langle m \rangle}(\mathbf{x}\_{2}), l\_{1}^{\langle m \rangle}(\mathbf{x}\_{3}), ..., l\_{1}^{\langle m \rangle}(\mathbf{x}\_{N}) \right]^{\top}, \\ \mathbf{u}^{\langle i \rangle} &=& \left[ u^{\langle i \rangle}(\mathbf{x}\_{1}), u^{\langle i \rangle}(\mathbf{x}\_{2}), u^{\langle i \rangle}(\mathbf{x}\_{3}), ..., u^{\langle i \rangle}(\mathbf{x}\_{N}) \right]^{\top}, \\ \mathbf{K}\_{1}^{\langle i \rangle} &=& \text{diag} \left( \kappa\_{1}^{\langle i \rangle}(\mathbf{x}\_{1}), \kappa\_{1}^{\langle i \rangle}(\mathbf{x}\_{2}), \kappa\_{1}^{\langle i \rangle}(\mathbf{x}\_{3}), ..., \kappa\_{1}^{\langle i \rangle}(\mathbf{x}\_{N}) \right). \end{array}$$

Then, we consider the second integral term with respect to space by letting it be *J m* <sup>2</sup> (*xk*) and using the idea of FIM-SCP (as described in Section 2.1) to approximate it. Then, we obtain

$$\mathrm{J}\_2^{\langle m\rangle}(\mathbf{x}\_k) := \int\_0^{\mathbf{x}\_k} \kappa\_2(\xi\_\nu^\mathbf{r}, t\_m) u(\xi\_\nu^\mathbf{r}, t\_m) d\xi = \int\_0^{\mathbf{x}\_k} \kappa\_2^{\langle m\rangle}(\xi) u^{\langle m\rangle}(\xi) d\xi \approx \sum\_{i=1}^N a\_{ki} \kappa\_2^{\langle m\rangle}(\mathbf{x}\_i) u^{\langle m\rangle}(\mathbf{x}\_i)$$

for each *xk* ∈ {*x*1, *x*2, *x*3, ..., *xN*}. The above equation can be written, in matrix form, as

$$\mathbf{J}\_2^{\langle m \rangle} = \mathbf{A} \mathbf{K}\_2^{\langle m \rangle} \mathbf{u}^{\langle m \rangle},\tag{13}$$

where **A** = **SS**¯ <sup>−</sup><sup>1</sup> is the shifted Chebyshev integration matrix defined in Section 2.1,

$$\begin{array}{rclclcl} \mathbf{J}\_{2}^{\langle m\rangle} &=& \left[J\_{2}^{\langle m\rangle}(\mathbf{x}\_{1}), J\_{2}^{\langle m\rangle}(\mathbf{x}\_{2}), J\_{2}^{\langle m\rangle}(\mathbf{x}\_{3}), \dots, J\_{2}^{\langle m\rangle}(\mathbf{x}\_{N})\right]^{\top}, \\ \mathbf{u}^{\langle m\rangle} &=& \left[u^{\langle m\rangle}(\mathbf{x}\_{1}), u^{\langle m\rangle}(\mathbf{x}\_{2}), u^{\langle m\rangle}(\mathbf{x}\_{3}), \dots, u^{\langle m\rangle}(\mathbf{x}\_{N})\right]^{\top}, \\ \mathbf{K}\_{2}^{\langle m\rangle} &=& \text{diag}\left(\kappa\_{2}^{\langle m\rangle}(\mathbf{x}\_{1}), \kappa\_{2}^{\langle m\rangle}(\mathbf{x}\_{2}), \kappa\_{2}^{\langle m\rangle}(\mathbf{x}\_{3}), \dots, \kappa\_{2}^{\langle m\rangle}(\mathbf{x}\_{N})\right). \end{array}$$

Then, we apply the FIM-SCP (described in Section 2.1) to eliminate all spatial derivatives from (11) by taking the *n*-layer integral on both sides of (11), to obtain the following equation at the shifted Chebyshev node *xk*, as defined in (2), as

$$\int\_0^{x\_k} \dots \int\_0^{\xi\_2} \left( \frac{u^{\langle m \rangle} - u^{\langle m - 1 \rangle}}{\pi} + \mathcal{L}u^{\langle m \rangle} \right) d\xi\_1 \dots d\xi\_n = \int\_0^{x\_k} \dots \int\_0^{\xi\_2} \left( J\_1^{\langle m \rangle} + J\_2^{\langle m \rangle} + F^{\langle m \rangle} \right) d\xi\_1 \dots d\xi\_n. \tag{14}$$

To simplify the *<sup>n</sup>*-layer integration of the spatial derivative terms of L*um*, by letting it be *<sup>Q</sup>m*(*xk*) and using the technique of integration by parts, we have

*Qm* (*xk*) := *xk* 0 ... *<sup>ξ</sup>*<sup>2</sup> 0 <sup>L</sup>*um* (*ξ*1) *dξ*1... *dξ<sup>n</sup>* = *xk* 0 ... *<sup>ξ</sup>*<sup>2</sup> 0 *n* ∑ *i*=0 *pi*(*ξ*1, *tm*) *<sup>d</sup><sup>i</sup> dx<sup>i</sup> um* (*ξ*1) *dξ*1... *dξ<sup>n</sup>* = *n* ∑ *i*=0 (−1)*<sup>i</sup> n i xk* 0 ... *<sup>η</sup>*<sup>2</sup> 0 *p* (*i*) *<sup>n</sup>* (*η*1, *tm*)*um* (*η*1)*dη*1... *dη<sup>i</sup>* + *xk* 0 *n*−1 ∑ *i*=0 (−1)*<sup>i</sup> n* − 1 *i <sup>ξ</sup><sup>n</sup>* 0 ... *<sup>η</sup>*<sup>2</sup> 0 *p* (*i*) *<sup>n</sup>*−1(*η*1, *tm*)*um* (*η*1)*dη*1... *dη<sup>i</sup> dξ<sup>n</sup>* + *xk* 0 *<sup>ξ</sup><sup>n</sup>* 0 *n*−2 ∑ *i*=0 (−1)*<sup>i</sup> n* − 2 *i <sup>ξ</sup>n*−<sup>1</sup> 0 ... *<sup>η</sup>*<sup>2</sup> 0 *p* (*i*) *<sup>n</sup>*−2(*η*1, *tm*)*um* (*η*1)*dη*1... *dη<sup>i</sup> dξn*−1*dξ<sup>n</sup>* . . . + *xk* 0 ... *<sup>ξ</sup>*<sup>2</sup> 0 *p*0(*ξ*1, *tm*)*um* (*ξ*1)*dξ*1... *dξ<sup>n</sup>* + *d*<sup>1</sup> *xn*−<sup>1</sup> *k* (*n* − 1)! + *d*<sup>2</sup> *xn*−<sup>2</sup> *k* (*n* − 2)! + *d*<sup>3</sup> *xn*−<sup>3</sup> *k* (*n* − 3)! + ... + *dn*,

where *d*1, *d*2, *d*3, ..., *dn* are the arbitrary constants which emerge from the process of integration by parts. Then, we substitute each *xk* ∈ {*x*1, *x*2, *x*3, ..., *xN*} into the above equation and utilize the idea of FIM-SCP. Thus, we can express it, in matrix form, by

$$\begin{split} \mathbf{Q}^{\langle m\rangle} &= \sum\_{i=0}^{n} (-1)^{i} \binom{n}{i} \mathbf{A}^{i} \mathbf{P}\_{n}^{(i)} \mathbf{u}^{\langle m\rangle} + \sum\_{i=0}^{n-1} (-1)^{i} \binom{n-1}{i} \mathbf{A}^{i+1} \mathbf{P}\_{n-1}^{\langle i\rangle} \mathbf{u}^{\langle m\rangle} \\ &+ \sum\_{i=0}^{n-2} (-1)^{i} \binom{n-2}{i} \mathbf{A}^{i+2} \mathbf{P}\_{n-2}^{\langle i\rangle} \mathbf{u}^{\langle m\rangle} + \dots + \mathbf{A}^{n} \mathbf{P}\_{0}^{\langle 0\rangle} \mathbf{u}^{\langle m\rangle} + \mathbf{X}\_{\text{n}} \mathbf{d} \\ &= \sum\_{j=0}^{n} \sum\_{i=0}^{n-j} (-1)^{i} \binom{n-j}{i} \mathbf{A}^{i+j} \mathbf{P}\_{n-j}^{\langle i\rangle} \mathbf{u}^{\langle m\rangle} + \mathbf{X}\_{\text{n}} \mathbf{d}, \end{split} \tag{15}$$

where **A** = **SS**¯ <sup>−</sup><sup>1</sup> is the shifted Chebyshev integration matrix, **d** = [*d*1, *d*2, *d*3, ..., *dN*] ,

$$\begin{array}{rcl} \mathbf{Q}^{\langle m\rangle} &=& \left[Q^{\langle m\rangle}(\mathbf{x}\_{1}), Q^{\langle m\rangle}(\mathbf{x}\_{2}), Q^{\langle m\rangle}(\mathbf{x}\_{3}), \dots, Q^{\langle m\rangle}(\mathbf{x}\_{N})\right]^{\top} \\ \mathbf{X}\_{\mathrm{il}} &=& \left[\mathbf{x}\_{\mathrm{n-1}}, \mathbf{x}\_{\mathrm{n-2}}, \mathbf{x}\_{\mathrm{n-3}}, \dots, \mathbf{x}\_{\mathrm{0}}\right] \text{ for each } \mathbf{x}\_{\mathrm{i}} = \frac{1}{\mathrm{l}!} \left[\mathbf{x}\_{1}^{\mathrm{i}}, \mathbf{x}\_{2}^{\mathrm{i}}, \mathbf{x}\_{3}^{\mathrm{i}}, \dots, \mathbf{x}\_{\mathrm{N}}^{\mathrm{i}}\right]^{\top} \\ \mathbf{P}^{\langle \mathbf{i}\rangle}\_{n-j} &=& \mathrm{diag}\left(p^{(\mathrm{i})}\_{n-j}(\mathbf{x}\_{1}, \mathbf{t}\_{m}), p^{(\mathrm{i})}\_{n-j}(\mathbf{x}\_{2}, \mathbf{t}\_{m}), p^{(\mathrm{i})}\_{n-j}(\mathbf{x}\_{3}, \mathbf{t}\_{m}), \dots, p^{(\mathrm{i})}\_{n-j}(\mathbf{x}\_{N}, \mathbf{t}\_{m})\right). \end{array}$$

Finally, we vary all points *xk* ∈ {*x*1, *x*2, *x*3, ..., *xN*} in (14) and rearrange them into matrix form by using the FIM-SCP with the derived matrix equations (12), (13), and (15); thus, we obtain

$$\frac{\mathbf{A}^n \mathbf{u}^{\langle m\rangle} - \mathbf{A}^n \mathbf{u}^{\langle m-1\rangle}}{\tau} + \mathbf{Q}^{\langle m\rangle} = \mathbf{A}^n \mathbf{J}\_1^{\langle m\rangle} + \mathbf{A}^n \mathbf{J}\_2^{\langle m\rangle} + \mathbf{A}^n \mathbf{F}^{\langle m\rangle}$$

or, factorizing the unknown solution **u***m* explicitly, as

$$\begin{cases} \mathbf{A}^{n} + \tau \sum\_{j=0}^{n} \sum\_{i=0}^{n-j} (-1)^{i} \binom{n-j}{i} \mathbf{A}^{i+j} \mathbf{P}\_{n-j}^{\langle i \rangle} - \frac{\tau^{2}}{2} \mathbf{A}^{n} \mathbf{K}\_{1}^{\langle m \rangle} - \tau \mathbf{A}^{n+1} \mathbf{K}\_{2}^{\langle m \rangle} \Big) \mathbf{u}^{\langle m \rangle} \\ \quad + \mathbf{X}\_{n} \mathbf{d} = \frac{\tau^{2}}{2} \mathbf{A}^{n} \mathbf{K}\_{1}^{\langle 0 \rangle} \mathbf{u}^{\langle 0 \rangle} + \tau^{2} \sum\_{i=1}^{m-1} \mathbf{A}^{n} \mathbf{K}\_{1}^{\langle i \rangle} \mathbf{u}^{\langle i \rangle} + \mathbf{A}^{n} \mathbf{u}^{\langle m-1 \rangle} + \tau \mathbf{A}^{n} \mathbf{F}^{\langle m \rangle}. \end{cases} \tag{16}$$

Next, consider the given boundary conditions (10) at the endpoints *b* ∈ {0, *L*}. We can convert them into matrix form by using the linear combination of shifted Chebyshev polynomial (4) in term of the *r*th-order derivative of *u* at the iteration time *tm* and using (3). Then, we have

$$\left. \frac{d^r}{dx^r} u^{\langle m \rangle} (x) \right|\_{x=b} = \sum\_{n=0}^{N-1} c\_n^{\langle m \rangle} \left. \frac{d^r}{dx^r} S\_n(x) \right|\_{x=b} = \psi\_\mathcal{I}(t\_m).$$

for all *r* ∈ {0, 1, 2, ..., *n* − 1}. We can express the above equation, in matrix form, as

$$
\begin{bmatrix} S\_0(b) & S\_1(b) & \cdots & S\_{N-1}(b) \\ S\_0'(b) & S\_1'(b) & \cdots & S\_{N-1}'(b) \\ \vdots & \vdots & \ddots & \vdots \\ S\_0^{(n-1)}(b) & S\_1^{(n-1)}(b) & \cdots & S\_{N-1}^{(n-1)}(b) \end{bmatrix} \begin{bmatrix} c\_0^{(m)} \\ c\_1^{(m)} \\ \vdots \\ c\_{N-1}^{(m)} \end{bmatrix} = \begin{bmatrix} \psi\_0(t\_m) \\ \psi\_1(t\_m) \\ \vdots \\ \psi\_{n-1}(t\_m) \end{bmatrix} \tag{17}
$$

which can be denoted by **Bc***m* = **Ψ***m* or **BS**−1**u***m* = **Ψ***m*. Finally, we can construct the system of *m*th iterative linear equations from (16) and (17), which has *N* + *n* unknowns containing **u***m* and **d**, as follows:

$$
\begin{bmatrix}
\mathbf{H}^{\langle m\rangle} & \mathbf{X}\_{\mathbf{n}} \\
\mathbf{B}\mathbf{S}^{-1} & \mathbf{0}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}^{\langle m\rangle} \\
\mathbf{d}
\end{bmatrix} = 
\begin{bmatrix}
\mathbf{E}\_{1}^{\langle m\rangle} \\
\mathbf{Y}^{\langle m\rangle}
\end{bmatrix}
\tag{18}
$$

where **H***m* is the coefficient matrix of **u***m* in (16) and **E***m* <sup>1</sup> is the right-hand side column vector of (16). Consequently, the solution **u***m* can be approximated by solving the system (18) starting from the given initial condition (9); that is, **u**0 = [*φ*(*x*1), *φ*(*x*2), *φ*(*x*3), ..., *φ*(*xN*)]. Note that, when we would like to find a numerical solution *u*(*x*, *t*) at any point *x* ∈ [0, *L*] for the terminal time *T*, we can calculate it by the following formula:

$$\mathfrak{u}(\mathfrak{x},T) = \sum\_{n=0}^{N-1} c\_n^{\langle m \rangle} S\_n(\mathfrak{x}) = \mathfrak{s}(\mathfrak{x}) \mathfrak{c}^{\langle m \rangle} = \mathfrak{s}(\mathfrak{x}) \mathbf{S}^{-1} \mathfrak{u}^{\langle m \rangle},$$

where **<sup>s</sup>**(*x*)=[*S*0(*x*), *<sup>S</sup>*1(*x*), *<sup>S</sup>*2(*x*), ..., *SN*−1(*x*)] and **<sup>u</sup>***m* is the final *<sup>m</sup>th* iterative solution of (18).
