(*Initial Parameters*)Needs["NDSolve`FEM`"];
Mobi = 1.0; lame = 0.01; noise = 0.02; conu0 = 0.63;
xmax = 1.0;
ymax = 1.0;
tmax = 1.0;

Ω = Rectangle[{0, 0}, {a, b}] /. {a -> 1, b -> 1};
RegionPlot[Ω, AspectRatio -> Automatic]
mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 1/1000, "MeshElementType" -> QuadElement];
mesh["Wireframe"]
n = Length[mesh["Coordinates"]]
u0 = ElementMeshInterpolation[{mesh}, conu0 + noise*(0.5 - RandomReal[{0, 1}, n])];
Plot3D[u0[x, y], {x, y} ∈ mesh]

op1 = D[u[t, x, y], t] - Laplacian[v[t, x, y], {x, y}] Mobi

op2 = v[t, x, y] - 200 u[t, x, y] (1 - 3 u[t, x, y] + 2 u[t, x, y]^2) +
      lame Laplacian[u[t, x, y], {x, y}]

{unn, vnn} =
  NDSolve[{op1 == 0, op2 == 0, u[0, x, y] == u0[x, y],
    v[0, x, y] == 0}, {u, v}, {t, 0, tmax}, {x, y} ∈ mesh];

(*InitialParameters*)
Needs["NDSolve`FEM`"];
Mobi = 1.0; lame = 0.01; noise = 0.02; conu0 = 0.63;
xmax = 1.0;
ymax = 1.0;
tmax = 1.0;
a = 1.;
b = 1.;

Ω = Rectangle[{0, 0}, {a, b}];
mesh = ToElementMesh[Ω,
   "MaxCellMeasure" -> {1 -> 0.005},
   "MeshElementType" -> QuadElement,
   "MeshOrder" -> 1
   ];

ClearAll[x, y, u];
vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {mesh}];
cdata = InitializePDECoefficients[vd, sd,
   "DiffusionCoefficients" -> {{-IdentityMatrix[2]}},
   "MassCoefficients" -> {{1}}
   ];
bcdata = InitializeBoundaryConditions[vd, sd, {{DirichletCondition[u[x, y] == 0., True]}}];
mdata = InitializePDEMethodData[vd, sd];

(*Discretization*)
dpde = DiscretizePDE[cdata, mdata, sd];
dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
{load, A, damping, M} = dpde["All"];
(*DeployBoundaryConditions[{load,A},dbc];*)
(*DeployBoundaryConditions[{load,M},dbc];*)

θ = 0.5;
τ = 0.000000001;
μ = Mobi;
λ = lame;
L = ArrayFlatten[{
    {M, τ μ θ A},
    {-λ A, M}
    }];
rhs[u_, v_] := Join[
   M.u - τ μ (1. - θ) A.v,
   M.(200. u (1. - 3. u + 2 u^2))
   ];
S = LinearSolve[L, Method -> "Pardiso"];

n = Length[mesh["Coordinates"]];
m = 10000;
f = x [Function] 100. ((1. - x^2)^2);
Df = x [Function] Evaluate[f'[x]];
rhs[u_, v_] := Join[M.u -  (μ τ (1. - θ)) A.v, M.Df[u]];
S = LinearSolve[L, Method -> "Pardiso"];

u0 = 2. RandomInteger[{0, 1}, n] - 1.;
ulist = ConstantArray[0., {m, n}];
ulist[[1]] = u = u0;

v0 = rhs[u0, 0. u0][[n + 1 ;; 2 n]];
v = v0;

Do[
  sol = S[rhs[u, v]];
  ulist[[k]] = u = sol[[1 ;; n]];
  v = sol[[n + 1 ;; 2 n]];
  , {k, 2, m}];

frames = Table[
   Image[
    Map[
     ColorData["ThermometerColors"],
     Partition[0.5 (Clip[ulist[[k]], {-1., 1.}] + 1.), Sqrt[n]],
     {2}
     ]
    ],
   {k, 1, m, 25}
   ];
Manipulate[
 frames[[k]],
 {k, 1, Length[frames], 1},
 TrackedSymbols :> {k}
 ]