This is the same approach that Tim describes in his solution, but implemented with the ToGradedMesh function added in version 13.0:

mesh = ToGradedMesh[
   Line[{{0}, {xmax}}], <|"Alignment" -> "Left", 
    "ElementCount" -> 100, "MinimalDistance" -> xmax/10000|>];
MeshRegion[mesh]

eqsH = {D[H[x, t], t] - dH D[H[x, t], x, x] == 
    NeumannValue[Ls Exp[a f Ea] P[x, t] - Ls Exp[-a f Ea] H[x, t], 
     x == 0], H[x, 0] == 1};
eqsP = {D[P[x, t], t] == 
    NeumannValue[-Ls Exp[a f Ea] P[x, t] + Ls Exp[-a f Ea] H[x, t], 
     x == 0], P[x, 0] == 1};

{Hfun, Pfun} = 
  NDSolveValue[{eqsH, eqsP}, {H, P}, x \[Element] mesh, {t, 0, tmax}, 
   Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}];

imgs = Plot[Hfun[x, #], x \[Element] mesh, 
     PlotRange -> {{0, 0.015}, {0.99999, 1.002}}] & /@ Subdivide[0, 2, 120];
ListAnimate@imgs

This is the same approach that Tim describes in his solution, but implemented with the ToGradedMesh function added in version 13.0:

mesh = ToGradedMesh[
   Line[{{0}, {xmax}}], <|"Alignment" -> "Left", 
    "ElementCount" -> 100, "MinimalDistance" -> xmax/10000|>];
MeshRegion[mesh]

eqsH = {D[H[x, t], t] - dH D[H[x, t], x, x] == 
    NeumannValue[Ls Exp[a f Ea] P[x, t] - Ls Exp[-a f Ea] H[x, t], 
     x == 0], H[x, 0] == 1};
eqsP = {D[P[x, t], t] == 
    NeumannValue[-Ls Exp[a f Ea] P[x, t] + Ls Exp[-a f Ea] H[x, t], 
     x == 0], P[x, 0] == 1};

{Hfun, Pfun} = 
  NDSolveValue[{eqsH, eqsP}, {H, P}, x \[Element] mesh, {t, 0, tmax}, 
   Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}];

imgs = Plot[Hfun[x, #], x \[Element] mesh, 
     PlotRange -> {{0, 0.015}, {0.99999, 1.002}}] & /@ Subdivide[0, 2, 120];
ListAnimate@imgs