f: array of numerical state at the interface between two cells
function crwenoL(n,u,f)
a = Array{Float64}(undef, n) # subdiagonal array of tridiagonal matrix
b = Array{Float64}(undef, n) # diagonal array of tridiagonal matrix
c = Array{Float64}(undef, n) # superdiagonal array of tridiagonal matrix
r = Array{Float64}(undef, n) # right hand~side
i=1
v1 = u[n-1]
v2 = u[n]
v3, v4, v5 = u[i], u[i+1], u[i+2]
a1,a2,a3,b1,b2,b3 = wcL(v1,v2,v3,v4,v5)
a[i], b[i], c[i] = a1, a2, a3
r[i] = b1*u[n] + b2*u[i] + b3*u[i+1]
i=2
v1 = u[n]
v2, v3, v4, v5 = u[i-1], u[i], u[i+1], u[i+2]
a1,a2,a3,b1,b2,b3 = wcL(v1,v2,v3,v4,v5)
a[i], b[i], c[i] = a1, a2, a3
r[i] = b1*u[i-1] + b2*u[i] + b3*u[i+1]
for i = 3:n-1
v1, v2, v3, v4, v5 = u[i-2], u[i-1], u[i], u[i+1], u[i+2]
a1,a2,a3,b1,b2,b3 = wcL(v1,v2,v3,v4,v5)
```

```
a[i], b[i], c[i] = a1, a2, a3
r[i] = b1*u[i-1] + b2*u[i] + b3*u[i+1]
end
i=n
v1, v2, v3, v4 = u[i-2], u[i-1], u[i], u[i+1]
v5 = u[2]
a1,a2,a3,b1,b2,b3 = wcL(v1,v2,v3,v4,v5)
a[i], b[i], c[i] = a1, a2, a3
r[i] = b1*u[i-1] + b2*u[i] + b3*u[i+1]
alpha = c[n]
beta = a[1]
ctdms(a,b,c,alpha,beta,r,f,1,n) # call cyclic Thomas algorithm function
end
```
For periodic boundary condition, we use the relation given by Equation (43) for ghost points. The tridiagonal system formed for the periodic boundary condition is not purely tridiagonal and it has the following form <sup>⎡</sup>

✝ ✆

$$
\begin{bmatrix} b\_1 & c\_1 & 0 & \cdots & \beta \\ a\_2 & b\_2 & c\_2 & \cdots & 0 \\ \vdots & \vdots & & \ddots & \vdots \\ 0 & \cdots & a\_{N-1} & b\_{N-1} & c\_{N-1} \\ \alpha & 0 & \cdots & a\_N & b\_N \end{bmatrix} \begin{bmatrix} u\_1 \\ u\_2 \\ \vdots \\ u\_{N-1} \\ u\_N \end{bmatrix} = \begin{bmatrix} d\_1 \\ d\_2 \\ \vdots \\ d\_{N-1} \\ d\_N \end{bmatrix},\tag{46}
$$

where *A***ˆ** is the tridiagonal system, *u* is the solution field and *d* is the right hand side of the tridiagonal system. We use Sherman-Morrison formula [17] to solve the periodic tridiagonal system given in Equation (46). The Sherman-Morrison formula can be used to solve a sparse linear system of the equation if there is a faster way of calculating *A*<sup>−</sup>1, where *A* is the modified matrix without *α* and *β*. We can compute the *A*−<sup>1</sup> using the regular Thomas algorithm explained in Section 2.3. We can write Equation (46) as given below

$$(A + p \otimes q)u = d,\tag{47}$$

where


where *<sup>γ</sup>* is arbitrary and the matrix *<sup>A</sup>* is the tridiagonal part of the matrix *<sup>A</sup>***ˆ**, and <sup>⊗</sup> represents the outer product of two vectors. Then the solution to Equation (46) can be obtained using below formulas

$$A \cdot y = d, \quad A \cdot z = p,\tag{49}$$

$$
\mu = y - \left[\frac{q \cdot y}{1 + q \cdot z}\right] z. \tag{50}
$$

The cyclic Thomas algorithm for solving the periodic tridiagonal system is given in Algorithm 2.

#### **Algorithm 2** Cyclic Thomas algorithm


Once the left and right side states (*uL*, *uR*) are available we compute the right hand side term of the Burgers equation using Listing 7. We use third-order Runge-Kutta numerical scheme for time integration. We solve the tridiagonal system formed with periodic boundary condition using cyclic Thomas algorithm and for Dirichlet boundary condition, we use regular Thomas algorithm. The implementation of cyclic Thomas algorithm in Julia is given in Listing 11. The implementation of regular Thomas algorithm is given in Listing 4.

**Listing 11.** Implementation of cyclic Thomas algorithm in Julia for periodic boundary condition domain. ✞ ☎

```