Mathematica code to solve numerically the Reaction - Diffusion PDEs
Leonardo Dagdug,  Ivan Pompa-García,  Jason Peña
From the book:
 Diffusion Under Confinement: A Journey Through Counterintuition

One-dimensional Gierer-Meinhardt model

Numerical solution of the coupled partial differential equations (PDEs) in  (13.55) and (13.56)  and its graphical representation.
The following code requires the values of the system length L,  reaction rate b, diffusion coefficient ratio d, and computation time tc,  in that specific order. Its output provides the numerical solution of u(x,t) and v(x,t), the concentration plot as a function of x, and a density plot of concentration u(x). For example, if we wish to see the results for a system with  L=30, b=0.35, d=30, and tc=400, we should write  PDEsSolution1D[30, 0.35, 30, 400].

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~IMPORTANT~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The entire code is confined within a single cell. It is imperative to exercise caution while executing the function, as it will yield the numerical solutions of the Partial Differential Equations (PDEs) and simultaneously generate both the concentration plots and the density plot
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~WARNING~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Be advised that the execution of the following command may necessitate a significant amount of time to complete, contingent upon the computational capabilities of your hardware.

(*Definition of the function PDEsSolution1D[L, b, d, tc] for the one-dimensional Gierer-Meinhardt reaction-diffusion model*)
PDEsSolution1D[L_ (*length of the system*),b_  (*reaction rate*), d_ (*difussivity ratio*),tc_ (*computation time*)]:=
(*Specification of the variables used within the Module*)
Module[{Dv=d, PDEs, initialconditions, x0, BCs, Eqs,Plotuv1D, Densityplotuv},

(******************************************NUMERICAL**SOLUTION**OF**THE**REACTION-DIFFUSION**SYSTEM**************************************)

(*First, we define the system of partial differential equations (PDEs) according to Eqs. (13.55) and (13.56), namely,*)
PDEs={(*Eq. (13.55)*)D[u[t,x],t]==Du D[u[t,x],{x,2}]+ (((u[t,x])^2)/(v[t,x]))-b u[t,x],
(*Eq.(13.56)*)D[v[t,x],t]== Dv D[v[t,x],{x,2}]+((u[t,x])^2)-v[t,x]};

(*The concentration gradients of the substances are represented as u[t,x] (activator) and v[t,x] (inhibitor).The variables inside the square brackets indicate the dependence of both functions u and v on time and space, while the command D[f[x1,x2,x3,...,xn],{xn,m}] gives the m-order partial derivative of function f with respect to xn. Thus, D[u[t,x],t] denotes partial derivative of u[t,x] with respect to time.*)
(*For the sake of simplicity, we set Du=1*)

Du=1;

(*Now, we set the initial and boundary conditions. On the one hand,a random value x0 is selected within the interval [0,L] to represent the initial position of species u, i.e.,*)

x0=RandomReal[{0,L}];                   

(*This selection is based on a Gaussian-like distribution centered at x0. On the other hand, we choose a constant value for v[0, x]. Therefore, the initial conditions are given by*)

initialconditions={ u[0,x]==1/b + Exp[-(((x-x0))^2)], v[0,x]==1/(b^2)};

(*When opting for periodic boundary conditions, we are lead to write*)

BCs={u[t,0]==u[t,L],v[t,0]==v[t,L]} ; 

(*Moving forward, we associate the initial and boundary conditions to the PDEs by using the Join command in Mathematica, this is*)

Eqs=Join[PDEs,initialconditions,BCs];

(*Finally, we solve numerically the system of PDEs over the given range [0, tc] for t and [0, L] for x. The solution is obtained by means of the NDSolve command and we name the solution solution1D*)

solution1D=NDSolve[Eqs (*System of PDEs to solve*),
{u,v} (*variables we wish to solve*),
{t,0,tc} (*computation time*),
{x,0,L} (*spatial domain*)];

(******************************************VISUALIZATION**OF**RESULTS**AND**DENSITY**PLOT**************************************)

(*The rest of the code is intended to show the position dependent plots of concentrations u(x) and v(x) in blue and orange,respectively, at different times. We incorporated controls to enable interactive time manipulation. Additionally, a density plot is displayed in the plane x-t.*)

(*Plot of the concentration gradients ~~the variable time is a dummy variable~~*)

Plotuv1D=
(*We use the Animate command to animate the plot*)
 Animate[
(*Evaluating functions u and v at the numerical solution solution1D and plotting its result*)Plot[Evaluate[{u[time,x], v[time,x]}/.solution1D],
(*Selecting spatial domain of the plot*)
 {x,0,L},
(*Detailed specification of legends*)
PlotLegends->{"u(x)", "v(x)"},
(*Frame Labels*)
Frame->True, FrameLabel->{"Position x", "Concentration gradient"} ,
(*Fixing plot range*)
PlotRange->{0, 12} ],
(*Time domain of the animation increasing in steps of tc/100*)
{time, 0,tc ,tc/100},                                    
(*Selecting manual reproduction of the animation*)
AnimationRunning->False];

(*Density plot of the concentration ~~the variable time is a dummy variable~~*)

Densityplotuv=  
(*Making a density plot of function u evaluated at the numerical solution solution1D*)
DensityPlot[Evaluate[u[time, x] /. solution1D],
(*Selecting spatial domain of the plot*)
{x, 0, L},
(*Time domain of the density plot*)
{time, 0, tc},
(*Plot label*)
PlotLabel->Style["Concentration u(x)"],
(*Detailed specification of frame labels*)
Frame->True,FrameLabel->{"Position x", "Time t"},
(*Definition of a color style*)
ColorFunction -> "SolarColors",
(*Using 100 initial sample points ~~Higher values give smoother plots~~*)
PlotPoints -> 100,
(*Size to display the plot*)
ImageSize->300] ;

(*Finally, we show all previous results in the form of a list, specifically the numerical solution, the concentration plot, and the density plot:*)

{solution1D, Plotuv1D, Densityplotuv}]

(*Subsequently, we undertake the solution of the system with L=80, b=0.35, d=30, and tc=400, using the latter function as PDEsSolution1D[80, 0.35, 30, 400]*)
PDEsSolution1D[100, 0.35, 30,400]

(*As we can see, the concentration gradient of both substances starts to oscillate along the x-axis and after a certain time both functions reach a steady-state*)


(*The reader may change the values of any of the parameters in order to see the role of the reaction rate b, critical length Lc, critical diffusion coefficient ratio dc, and time computation (Subscript[t, c]). For instance, by using Eq. (13.74), we see that the critical lenght is about Lc= 12.0, meaning that below that value, the concentration remain approximately constant and no patterns emerge in the density plot. Therefore, as the domain size increases,the substance concentration exhibits steeper transitions*)