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

Two-dimensional Gierer-Meinhardt model

Numerical solution of the coupled partial differential equations (PDEs) in  (13.73) and (13.74)  and its graphical representation.
The following code requires the values of the system length Lx and Ly,  reaction rate b, diffusion coefficient ratio d, and computation time tc,  in that specific order. Its output provides the numerical solution of u(x,y,t) and v(x,y,t), the concentration plot as a function of x and y, and a density plot of concentrations. For example, if we wish to see the results for Lx=80, Ly=10, b=0.45, d=30, and tc=1000, we should write  PDEsSolution2D[80, 10, 0.45, 30, 1000].

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~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  the density plot animation.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~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 PDEsSolution2D[Lx,Ly,b, d, tc] for the two-dimensional Gierer-Meinhardt reaction-diffusion model*)

PDEsSolution2D[Lx_ (*length along the x-axis*),Ly_ (*length along the y-axis*),
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, Densityplotuv},

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

(*First, we define the system of partial differential equations (PDEs) according to Eqs. (13.73) and (13.74), namely,*)

PDEs={(*Eq.(13.73)*)D[u[t,x,y], t]==Du (D[u[t,x,y],{x,2}]+D[u[t,x,y],{y,2}])+(((u[t,x, y])^2)/(v[t,x, y]))-b u[t,x, y], 
(*Eq.(13.74)*)D[v[t,x,y], t]==Dv (D[v[t,x,y],{x,2}]+D[v[t,x,y],{y,2}])+((u[t,x, y])^2)-v[t,x, y]};

(*The concentration gradients of the substances are represented as u[t,x,y] (activator) and v[t,x,y] (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,...,xn],{xn,m}] gives the m-order partial derivative of function f with respect to xn. Thus, D[u[t,x,y],t] denotes partial derivative of u[t,x,y] 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 values x0 is selected within the interval [0, Lx] to represent the initial position of species u . This selection is based on a Gaussian-like distribution centered at (x0, y)., i.e.,*)

x0=RandomReal[{0, Lx}];

 (*On the other hand, we choose a constant value for v[0,x,y]. Therefore, the initial conditions are given by*)

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

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

BCs={u[t,0,y]==u[t,Lx,y],v[t,x,0]==v[t, x, Ly]}; 

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

 Eqs=Join[PDEs,initialconditions,BCs];

(*Finally, we solve numerically the system of PDEs in the range [0, tc] for t, [0, Lx] for x,and [0, Ly] for y. The solution is obtained by meansthe NDSolve comman and we name the solution solution2D*)

solution2D=NDSolve[Eqs (*System of PDEs to solve*),{u,v},
{t,0,tc} (*computation time*),
{y,0,Ly}(*spatial domain along the y-axis*),
{x,0,Lx} (*spatial domain along the x-axis*)];

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

(*The rest of the code is intended to show the density plots of concentrations u(x, y) and v(x, y) with incorporated controls to enable interactive time manipulation.*)

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

Densityplotuv=   
(*We use the ListAnimate command to animate results of the Table below*)
ListAnimate[
(*Generating a list of the DensityPlot when time goes from 0 to tc*)
Table[ 
(*Making a density plot of function v evaluated at the numerical solution solution2D*)
 DensityPlot[Evaluate[v[time, x, y] /. solution2D],
(*Selecting spatial domain of the plot*)
 {x, 0, Lx},{y, 0, Ly},
(*Definition of a color style*)
               ColorFunction -> "SolarColors",
(*Plot label*)
PlotLabel->Style["Concentración de u(x,t)"],
(*Frame label*)
FrameLabel->{"Position x","Position y"},
(*Using 100 initial sample points ~~Higher values give smoother plots~~*)
 PlotPoints -> 100,
(*Automatic set up of ratio of height to width of the plot*)
 AspectRatio->Automatic],
(*Time domain of the animation increasing in steps of tc/100*)
 {time, 0, tc, tc/100}],
(*Selecting manual reproduction of the animation*)
 AnimationRunning->False];


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

{solution2D,  Densityplotuv}]

(*Subsequently, we undertake the solution of the system with Lx=80, Ly=10, b=0.45, d=30, and tc=1000, using the latter function as PDEsSolution2D[80, 10, 0.45, 30, 1000]*)

PDEsSolution2D[80  (*length along the x-axis*),80  (*length along the y-axis*), 0.45 (*reaction rate*), 30 (*difussivity ratio*),
2000 (*computation time*)]


(*The reader may change the values of any of the parameters in order to see the role of the reaction rate b, domain size, critical diffusion coefficient ratio dc, and time computation tc. For instance, when using PDEsSolution2D[80, 10, 0.35, 30, 400], the resultant patterns show stripes, whereas upon employing a different configuration, say a square system with PDEsSolution2D[80, 80, 0.45, 30, 1000], dots are observed.  As in the one-dimensional case, we conclude that the larger the diffusion reaction domain, the more complex the spatial patterns.*)