un: matrix for storing solution at every time step
for k = 2:nt+1
for i = 2:nx
un[k,i] = un[k-1,i] + beta*(un[k-1,i+1] - 2.0*un[k-1,i] + un[k-1,i-1])
end
un[k,1] = 0.0 # boundary condition at x = -1
un[k,nx+1] = 0.0 # boundary condition atx=1
end
✝ ✆
```
The comparison of the exact and numerical solution using the FTCS scheme and absolute error plot are shown in Figure 2. We see that we get a very good agreement between the exact and numerical solution due to very small grid size in both space and temporal direction. If we use a larger grid size, then we will see deprecation in accuracy.

**Figure 2.** Comparison of exact solution and numerical solution computed using FTCS numerical scheme. The numerical solution is computed using Δ*t* = 0.0025 and Δ*x* = 0.025.

#### *2.2. Runge-Kutta Numerical Scheme*

We saw in Section 2.1 that the accuracy of the numerical scheme depends upon the number of terms included from the Taylor series expansion. One of the ways to increase accuracy is to include more terms. The additional term involves partial derivative of *f*(*x*, *t*) which provides additional information on *f* at *t* = *tn*. For the time stepping, we only have information about the data at one or more time steps and the analytical form of *f* is not known between two time steps. Runge-Kutta methods tries to improve the accuracy of temporal term by evaluating *f* at intermediate points between *tn* and *tn*+1. The additional steps lead to an increase in computational time, but the temporal accuracy is increased.

We use a third-order accurate Runge-Kutta scheme [16] for the discretization of the temporal term in the heat equation. We use the same second-order central difference scheme for the spatial term. The truncation error of this numerical approximation of heat equation is O(Δ*t* 3, Δ*x*2). We move from time step *tn* to *tn*+<sup>1</sup> using three steps. The time integration of the heat equation using third-order Runge-Kutta scheme is given below:

$$u\_i^{(1)} = u\_i^{(n)} + \frac{a\Delta t}{\Delta x^2} (u\_{i+1}^{(n)} - 2u\_i^{(n)} + u\_{i-1}^{(n)}),\tag{12}$$

$$u\_i^{(2)} = \frac{3}{4}u\_i^{(n)} + \frac{1}{4}u\_i^{(1)} + \frac{a\Delta t}{4\Delta x^2} (u\_{i+1}^{(n)} - 2u\_i^{(n)} + u\_{i-1}^{(n)}),\tag{13}$$

$$u\_i^{(n+1)} = \frac{1}{3}u\_i^{(n)} + \frac{2}{3}u\_i^{(2)} + \frac{2a\Delta t}{3\Delta x^2} (u\_{i+1}^{(n)} - 2u\_i^{(n)} + u\_{i-1}^{(n)}).\tag{14}$$

We can also construct the higher-order Runge-Kutta numerical scheme using more intermediate steps. The implementation of the third-order Runge-Kutta numerical scheme in Julia is given in Listing 2. The comparison of the exact and numerical solutions using the Runge-Kutta scheme and absolute error plot is shown in Figure 3. The discretization error is less for Runge-Kutta numerical scheme compared to the FTCS numerical scheme due to a higher-order approximation of temporal term.

```
Listing 2. Implementation of third-order Runge-Kutta numerical scheme for heat equation in Julia. ✞ ☎
```

```