Here is a link to a SageMathCell implementation for the discrete logistic model reaction-diffusion system. I also give the code here:

import scipy;
from scipy import integrate
from sage.plot.colors import rainbow

pts = 20 # number of spatial points
steps = 200 # time steps
deltat = 0.1 # size of time step
t = srange(0, steps*deltat, deltat) # list of time values

# Neumann matrix
NeuM = -2*identity_matrix(pts)
for k in [0..pts-2]:
    NeuM[k+1,k]=1
    NeuM[k, k+1]=1
NeuM[0,0]=-1
NeuM[pts-1, pts-1] = -1

# function to square a list
def square(x):
    for k in [0..len(x)-1]:
        x[k]=x[k]*x[k]
    return(x)

# right hand side of our system
mu = 0.5 # diffusion constant
def F(x, t=0):
    return(mu*NeuM*vector(x) + x - square(x))

# initial condition all zero
u0 = [0] 
for k in [1..pts-1]:
    u0 = u0 + [0]

# except for a small spike
u0[1]=.01
    
# evolve using odeint (switch to RK4 at some point?)
u = integrate.odeint(F, u0, t)

#plot each location in corresponding color
colorz = rainbow(pts)
mainplot = line(zip(t,u.transpose()[0]), color = colorz[0])
for k in[1..pts-1]:
    mainplot += line(zip(t,u.transpose()[k]), color = colorz[k])
mainplot += line(zip(t,u.transpose()[1]), color = "black")


# plot spatial profile at each time
colorz = rainbow(steps)
profile = list_plot(u[0], plotjoined=true, color = colorz[0])
for k in [1..steps-1]:
    profile +=list_plot(u[k], plotjoined=true, color = colorz[k])
profile +=list_plot(u[0], plotjoined=true, color = "black")

# display plots for each location
mainplot.show(figsize=[4,3], ymax = 1.2)
# output profile time series as pdf
profile.save(filename="profile.pdf")