Simulating a partial differential equation, uniform vegetation and bare soil system (2 species)

Question

I want to simulate a system of two plants ,woody and herbaceous species; with starting condition that half of the grid (square grid) is filled with woody specie, the herbaceous specie exist on the line between the woody specie and the bare soil.

The system model is the following:

$\frac{\partial B_1}{\partial t}= \Lambda_1 WB_1(1-B_1/K_1)(1+E_1B_1)^2-M_1B_1+D_{B_1}\nabla^2B_1 ,$

$\frac{\partial B_2}{\partial t}= \Lambda_2 WB_2(1-B_2/K_2)(1+E_2B_2)^2-M_2B_2+D_{B_2}\nabla^2B_2 ,$

$\frac{\partial W}{\partial t} = I H - \frac{N W}{1+RB_1/K_1} -\Gamma_1WB_1(1+E_1B_1)^2-\Gamma_2WB_2(1+E_2B_2)^2+D_W\nabla^2W,$

$\frac{\partial H}{\partial t}=P-IH +D_H\nabla^2(H^2)$

where

$I=A\frac{Y_1B_1+Y_2B_2+qf}{Y_1B_1+Y_2B_2+q}.$

subject to Newman Boundary Conditions and Initial Conditions i mentioned above, I wrote the following code, but with no success.

L = 1;(*length of square*)
T = 8;(*Time integration*)
pts = 100;
f = 0.2; (*The infiltration contrast between vegetated and bare-soil \
areas*)
A = 40.0; (*The maximal infiltration rate, obtained in densely \
vegetated areas*)
n = 30.0; (*the soil-water evaporation rate*)
R = 30.0; (*the shading effect*)
E1 = 2.0; (*root to shoot ratio of B1*)
E2 = 0.0;(*root to shoot ratio of B2*)
M1 = 0.75; (*the mortality rate of B1*)
M2 = 7.5; (*the mortality rate of B2*)
K1 = 5.0; (*the maximum standing biomass of B1*)
K2 = 0.5; (*the maximum standing biomass of B2*)
\[CapitalLambda]1 = 0.05; (*growth rate of B1*)
\[CapitalLambda]2 = 1.0; (*growth rate of B2*)
\[CapitalGamma]1 = 10.0; (*the water uptake rate of B1*)
\[CapitalGamma]2 = 15.0; (*the water uptake rate of B2*)
q = 5.0; (*reference biomass value beyond which he effect of \
vegetation on infiltration becomes significant*)
Y1 = 1.0; (*relative contribution of the woody specie (B1)*)
Y2 = 100; (*relative contribution of the herbaceous specie (B2)*)
DB1 = 0.1; (*the biomass expansion rate of B1*)
DB2 = 0.05; (*the biomass expansion rate of B2*)
DW = 1.0; (*the soil-water diffusion*)
DH = 1.0; (*the rate of surface-water spread*)
i = A (Y1 B1[t, x, y] + Y2 B2[t, x, y] + q f)/(
  Y1 B1[t, x, y] + Y2 B2[t, x, y] + q );
P = 2;
(*system of nonlinear PDE*)
pde = {D[B1[t, x, y], 
    t] == \[CapitalLambda]1 W[t, x, y] B1[t, x, 
      y] (1 - B1[t, x, y]/K1) (1 + E1 B1[t, x, y])^2 - 
    M1 B1[t, x, y] + 
    DB1 (D[B1[t, x, y], x, x] + D[B1[t, x, y], y, y]), 
  D[B2[t, x, y], 
    t] == \[CapitalLambda]2 W[t, x, y] B2[t, x, 
      y] (1 - B2[t, x, y]/K2) (1 + E2 B2[t, x, y])^2 - 
    M2 B2[t, x, y] + 
    DB2 (D[B2[t, x, y], x, x] + D[B2[t, x, y], y, y]),
  D[W[t, x, y], t] == 
   i H - n W[t, x, y]/(
     1 + R B1[t, x, y]/K1 ) - \[CapitalGamma]1 W[t, x, y] B1[t, x, 
      y] (1 + E1 B1[t, x, y])^2 - \[CapitalGamma]2 W[t, x, y] B2[t, 
      x, y] (1 + E2 B2[t, x, y])^2 + 
    DW (D[W[t, x, y], x, x] + D[W[t, x, y], y, y]),
  D[H[t, x, y], t] == 
   P - i H[t, x, y] + 
    DH (D[H[t, x, y], x, x]^2 + D[H[t, x, y], y, y]^2)};
(*Newman boundary condition*)
bc = {(D[B1[t, x, y], x] /. x -> -L) == 
   0, (D[B1[t, x, y], x] /. x -> L) == 
   0, (D[B1[t, x, y], y] /. y -> -L) == 
   0, (D[B1[t, x, y], y] /. y -> L) == 
   0, (D[B2[t, x, y], x] /. x -> -L) == 
   0, (D[B2[t, x, y], x] /. x -> L) == 
   0, (D[B2[t, x, y], y] /. y -> -L) == 
   0, (D[B2[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,
  (D[H[t, x, y], x] /. x -> -L) == 0, (D[H[t, x, y], x] /. x -> L) ==
    0, (D[H[t, x, y], y] /. y -> -L) == 
   0, (D[H[t, x, y], y] /. y -> L) == 0};

ic = {B1[0, x, y] == 
   Interpolation[
     Flatten[Table[{x, y, 1}, {x, -L, L, 2/pts}, {y, -L, 0, 2/pts}], 
      1]][x, y], 
  B1[0, x, y] == 
   Interpolation[
     Flatten[Table[{x, y, 0}, {x, -L, L, 2/pts}, {y, 0, L, 2/pts}], 
      1]][x, y],
  B2[0, x, y] == 
   Interpolation[
     Flatten[Table[{x, y, 1}, {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],
  H[0, x, y] == 
   Interpolation[
     Flatten[Table[{x, y, RandomReal[]}, {x, -L, L, 2/pts}, {y, -L, 
        L, 2/pts}], 1]][x, y]};

eqns = Flatten@{pde, bc, ic};

sol = NDSolve[
  eqns, {B1[t, x, y], B2[t, x, y], W[t, x, y], H[t, x, y]}, {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 + 3*t1;
 DensityPlot[B1[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]

Can someone help me? Thanks in advance.