When plotting the numerical solution at step 4, initial conditions of the equation was taken from analytical solution (by using un = u.copy() at the beginning of the FD loop, array u is the analytical solution of function at t=0 in our code). I think we need to use initial conditions below in the picture.
