Solving a Reaction-diffusion system

I'm just a beginner in Mathematica software. I have version 12.

I am looking for a program that will solve a prey-predator system with diffusion. My problem is that I tried a previously posted program, but the code is not compatible with Mathematica 12. Another problem is that I need to plot the solution as an animation to show the periodic motion. The program is next:

 pts = 100;
t0 = 0;
t1 = 2;
tmax = 3000;
(*length of square*)  L = 1;(*Time integration*) T = 8;(* Diffusion \
parameter for the prey*)d1 = 0.00028;(* Diffusion parameter for the \
predator*)d2 = 0.00028;(* Fertility parameter for the prey*)a = 1;(* \
Mortality parameter of the prey in the presence of predator*)b = 1;(* \
Fertility parameter of the predator*) c = 1;(* Fertility parameter of \
the predator in the presence of the prey*) d = 1;
(* Mortality parameter of the predator*) e = 1;

(*system of nonlinear PDE*)
pde = {D[u[t, x, y], t] ==
d1 (D[u[t, x, y], x, x] + D[u[t, x, y], y, y]) + a u[t, x, y] -
b u[t, x, y]*w[t, x, y],
D[w[t, x, y], t] ==
d2 (D[w[t, x, y], x, x] + D[w[t, x, y], y, y]) + c w[t, x, y] +
d u[t, x, y]*w[t, x, y] - e w[t, x, y]};

(*Newman boundary condition*)

bc = {(D[u[t, x, y], x] /. x -> -L) ==
0, (D[u[t, x, y], x] /. x -> L) ==
0, (D[u[t, x, y], y] /. y -> -L) ==
0, (D[u[t, x, y], y] /. y -> L) ==
0, (D[w[t, x, y], x] /. x -> -L) ==
0, (D[w[t, x, y], x] /. x -> L) ==
0, (D[w[t, x, y], y] /. y -> -L) ==
0, (D[w[t, x, y], y] /. y -> L) == 0};
(*initial condition*)

ic = {u[0, x, y] ==
Interpolation[
Flatten[Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L,
L, 2/pts}], 1]][x, y],
w[0, x, y] ==
Interpolation[
Flatten[Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L,
L, 2/pts}], 1]][x, y]};
(*ic={u[0,x,y]\[Equal]0,w[0,x,y]\[Equal]0};*)

eqns = Flatten@{pde, bc, ic};
sol = NDSolve[eqns, {u, w}, {t, 0, T}, {x, -L, L}, {y, -L, L},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> pts, "MaxPoints" -> pts},
Method -> {"IDA", "ImplicitSolver" -> {"GMRES"}}}];

GraphicsGrid[Table[time = t0 + 2*t1;
DensityPlot[u[time, x, y] /. sol, {x, -L, L}, {y, -L, L},
ColorFunction -> "SunsetColors",
PlotLabel -> "t=" <> ToString[time], Ticks -> False], {t1, 0,
2}, {t0, 0, 2}], ImageSize -> 600]
(Animate[DensityPlot[
Evaluate[u[x, y, time] /. sol], {x, 0, L}, {y, 0, L},
FrameLabel -> {"x", "y"}, PlotPoints -> pts,
ColorFunctionScaling -> False], {time, 0, 3000}]


It appears the following warnings: NDSolve: Warning: boundary and initial conditions are inconsistent NDSolve: Warning: estimated initial error on the specified spatial grid in the direction of independent variable x exceeds prescribed error tolerance NDSolve: "Further output of

 \!$$\*StyleBox[RowBox[{\"InterpolatingFunction\", \ \"::\", \"dmval\"}], \"MessageName\"]$$ will be suppressed during \


this calculation."

Also the plot anime doesn't out.

If someone will help me, I will be very grateful.

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, [by clicking the checkmark sign](tinyurl.com/4srwe26 – Dunlop Jan 7 at 19:46
• Have you seen this? – Henrik Schumacher Jan 7 at 19:47
• Perhaps to get more assistance you can show the code (plus the link to the code) that you used. Furthermore there are quite some examples of similar codes that you can find if you search for Reaction-Diffusion equations – Dunlop Jan 7 at 19:48
• Does this answer your question? Simulating a partial differential equation - reaction-diffusion systems and Turing patterns – Dunlop Jan 7 at 19:50
• Kamal, the example form (178559) does run in version 12. It just throws warnings because several parameters (for example the boundary conditions and the time step size) have been chosen suboptimally. – Henrik Schumacher Jan 8 at 9:37

The Helmholtz equation is a reaction-diffusion equation. There are many examples in the documentation.

Maybe something like this:

region = RegionDifference[Rectangle[{0, 0}, {5, 10}], Disk[{5, 5}, 3]];
op = -Laplacian[u[x, y], {x, y}] - u[x, y];
bcs = {DirichletCondition[u[x, y] == 0, x == 0 && 8 <= y <= 10],
DirichletCondition[u[x, y] == 100, (x - 5)^2 + (y - 5)^2 == 3^2]};
ufun = NDSolveValue[{op == 0, bcs}, u, Element[{x, y}, region]];
ContourPlot[ufun[x, y], Element[{x, y}, region],
ColorFunction -> "Temperature", AspectRatio -> Automatic]