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")