Can you please revise the python code import numpy as np import matplotlib.pyplot as plt N = 20 # number of points to discretize L = 1.0 X = np.linspace(0, L, N) # position along the rod h = L / (N - 1) # discretization spacing C0t = 0.1 # concentration at x = 0 D = 0.02 tfinal = 50.0 Ntsteps = 1000 dt = tfinal / (Ntsteps - 1) t = np.linspace(0, tfinal, Ntsteps) alpha = D * dt / h**2 print (alpha) C_xt = [] # container for all the time steps # initial condition at t = 0 C = np.zeros(X.shape) C[0] = C0t C_xt += [C] for j in range(1, Ntsteps): N = np.zeros(C.shape) N[0] = C0t N[1:-1] = alpha*C[2:] + (1 - 2 * alpha) * C[1:-1] + alpha * C[0:-2] N[-1] = N[-2] # derivative boundary condition flux = 0 C[:] = N C_xt += [N] # plot selective solutions if j in [1,2,5,10,20,50,100,200,500]: plt.plot(X, N, label='t={0:1.2f}'.format(t[j])) plt.xlabel('Position in rod') plt.ylabel('Concentration') plt.title('Concentration at different times') plt.legend(loc='best') plt.savefig('transient-diffusion-temporal-dependence.png') C_xt = np.array(C_xt) plt.figure() plt.plot(t, C_xt[:,5], label='x={0:1.2f}'.format(X[5])) plt.plot(t, C_xt[:,10], label='x={0:1.2f}'.format(X[10])) plt.plot(t, C_xt[:,15], label='x={0:1.2f}'.format(X[15])) plt.plot(t, C_xt[:,19], label='x={0:1.2f}'.format(X[19])) plt.legend(loc='best') plt.xlabel('Time') plt.ylabel('Concentration') plt.savefig('transient-diffusion-position-dependence.png') plt.show()

Can you please revise the python code 

import numpy as np

import matplotlib.pyplot as plt

N = 20  # number of points to discretize

L = 1.0

X = np.linspace(0, L, N) # position along the rod

h = L / (N - 1)          # discretization spacing

C0t = 0.1  # concentration at x = 0

D = 0.02

tfinal = 50.0

Ntsteps = 1000

dt = tfinal / (Ntsteps - 1)

t = np.linspace(0, tfinal, Ntsteps)

alpha = D * dt / h**2

print (alpha)

C_xt = [] # container for all the time steps

# initial condition at t = 0

C = np.zeros(X.shape)

C[0] = C0t

C_xt += [C]

for j in range(1, Ntsteps):

    N = np.zeros(C.shape)

    N[0] =  C0t

    N[1:-1] = alpha*C[2:] + (1 - 2 * alpha) * C[1:-1] + alpha * C[0:-2]

    N[-1] = N[-2]  # derivative boundary condition flux = 0

    C[:] = N

    C_xt += [N]

    # plot selective solutions

    if j in [1,2,5,10,20,50,100,200,500]:

        plt.plot(X, N, label='t={0:1.2f}'.format(t[j]))

plt.xlabel('Position in rod')

plt.ylabel('Concentration')

plt.title('Concentration at different times')

plt.legend(loc='best')

plt.savefig('transient-diffusion-temporal-dependence.png')

C_xt = np.array(C_xt)

plt.figure()

plt.plot(t, C_xt[:,5], label='x={0:1.2f}'.format(X[5]))

plt.plot(t, C_xt[:,10], label='x={0:1.2f}'.format(X[10]))

plt.plot(t, C_xt[:,15], label='x={0:1.2f}'.format(X[15]))

plt.plot(t, C_xt[:,19], label='x={0:1.2f}'.format(X[19]))

plt.legend(loc='best')

plt.xlabel('Time')

plt.ylabel('Concentration')

plt.savefig('transient-diffusion-position-dependence.png')

plt.show()