un: matrix for storing solution at every time step
s, e = 1, nx
for k = 2:nt+1
i=1 # left boundary
a[i], b[i], c[i], d[i] = 0.0, 1.0, 0.0, 0.0
for i = 2:nx # interior points
a[i], b[i], c[i] = -beta, 1.0+2.0*beta, -beta
d[i] = beta*un[k-1,i+1] + (1.0-2.0*beta)*un[k-1,i] + beta*un[k-1,i-1]
end
i = nx+1 # right boundary
a[i], b[i], c[i], d[i] = 0.0, 1.0, 0.0, 0.0
tdms(a,b,c,r,p,s,e) # call Thomas algorithm function to calculate p
un[k,:] = p
end
```
**Figure 4.** Comparison of exact solution and numerical solution computed using Crank-Nicolson numerical scheme. The numerical solution is computed using Δ*t* = 0.0025 and Δ*x* = 0.025.

**Listing 4.** Implementation of Thomas algorithm to solve tridiagonal system of equations in Julia. ✞ ☎

```