Specifying Prescribed Mesh to Solve for a System Consisting of Two PDEs and an ODE

Question

I am solving the following equations: $$ \frac{\delta C_{1\text{f}}^*}{\delta t^*} = \frac{\overline{t}}{\overline{x}^2} \frac{\delta}{\delta x^*} \left\{ \mathcal{D} \left[ \frac{\delta C^*_{1\text{f}}}{\delta x^*} + \beta z_1 e C^*_{1\text{f}} \frac{\delta \psi}{\delta x^*} \right] \right\} $$ $$ \frac{\delta C_{2\text{f}}^*}{\delta t^*} = \frac{\overline{t}}{\overline{x}^2} \frac{\delta}{\delta x^*} \left\{ \mathcal{D} \left[ \frac{\delta C^*_{2\text{f}}}{\delta x^*} + \beta z_1 e C^*_{2\text{f}} \frac{\delta \psi}{\delta x^*} \right] \right\} - \frac{\delta C_{2\text{b}}^*}{\delta t^*} $$ $$ \frac{\delta}{\delta x^*} \left( \frac{\epsilon \epsilon_0}{\overline{x}^2 e N_A} \frac{\delta \psi}{\delta x^*}\right) = - ( -C_F +z_2(C_{2\text{b}}+C_{2\text{f}}) +z_1 C_1) $$ $$ C_{2\text{b}}^* = C_{F}^* \left( 1 +\exp{\beta z_2 \psi} \right)^{-1} $$ With initial conditions: $$ \begin{equation} C_{1\text{f}}^*(x^*,0) = \begin{cases} \overline{C} & x^*<3.5,\\ 0 & x^*>3.5 \end{cases} \end{equation} $$ $$ \begin{equation} C_{2\text{f}}^*(x^*,0) = \begin{cases} -\frac{z_2}{z_1} \overline{C} & x^*<3.5,\\ 0 & x^*>3.5 \end{cases} \end{equation} $$ And boundary conditions: $$ \begin{equation} C_{1\text{f}}^*(0,t^*) = -\frac{z_2}{z_1} \overline{C} \end{equation} $$ $$ \begin{equation} \frac{\delta C_{1\text{f}}^* (4,t^*)}{\delta x^*} = 0 \end{equation} $$ $$ \begin{equation} C_{2\text{f}}^*(0,t^*) = \overline{C} \end{equation} $$ $$ \begin{equation} \frac{\delta C_{2\text{f}}^* (4,t^*)}{\delta x^*} = 0 \end{equation} $$ $$ \begin{equation} \psi(0,t^*) = 0 \end{equation} $$ $$ \begin{equation} \frac{\delta \psi(4,t^*)}{\delta x^*} = 0 \end{equation} $$ My initial codes are as follows:

ClearAll["Global`*"]
(*Parameters*)
kB = 1.380649*10^-23; 
absTemp = 298.15;
beta = 1/(kB*absTemp);
eps0 = 8.8541878128*10^(-12);
eps = 78.6;
e = 1.60217663*10^(-19);
nA = 6.02214076*10^(23);
dfsvt = 1.33*10^(-9);
xRef = 10*10^(-6);
tRef = 10*10^(-2);
xStarMem0 = 3;
xStarMemMid = 3.5;
cSaltBar = 1500;(*mol/m3*)
cFVal = 3000;(*mol/m3*)
cF[xStar_] = cFVal/2*(Tanh[15*(xStar - xStarMem0)] + 1); 
z1 = -1;
z2 = 1;

(*BC and IC*)
c1fIC = c1f[xStar, 0] == -z2/z1*
    cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1); 
c1fDirichBC = c1f[0, tStar] == -z2/z1*cSaltBar;
c1fNuemBC = NeumannValue[0, xStar == xStarMemMid];
c2fIC = c2f[xStar, 0] == 
   cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1); 
c2fDirichBC = c2f[0, tStar] == cSaltBar;
c2fNuemBC = NeumannValue[0, xStar == xStarMemMid];
phiDirichBC = phi[0, tStar] == 0;
phiNuemBC = NeumannValue[0, xStar == xStarMemMid];

(*Equations*)
c2b[xStar_, tStar_] = 
  cF[xStar]*(1 + Exp[beta*z2*e*phi[xStar, tStar]])^(-1);
c1Trnsprt = 
  D[c1f[xStar, tStar], tStar] == 
   tRef/(xRef^2)*
     D[dfsvt*(D[c1f[xStar, tStar], xStar] + 
         beta*z1*e*c1f[xStar, tStar]*D[phi[xStar, tStar], xStar]), 
      xStar] + c1fNuemBC;
c2Trnsprt = 
  D[c2f[xStar, tStar], tStar] == 
   tRef/(xRef^2)*
     D[dfsvt*(D[c2f[xStar, tStar], xStar] + 
         beta*z2*e*c2f[xStar, tStar]*D[phi[xStar, tStar], xStar]), 
      xStar] - D[c2b[xStar, tStar], tStar] + c2fNuemBC;
elctrcTrnsprt = 
  D[eps*eps0/(xRef^2*e*nA)*D[phi[xStar, tStar], xStar], 
    xStar] == -(-cF[xStar] + 
       z2*(c2f[xStar, tStar] + c2b[xStar, tStar]) + 
       z1*c1f[xStar, tStar]) + phiNuemBC;

(*Input Parameter Plotting*)
Plot[cF[x], {x, 0, xStarMemMid}, PlotRange -> All, 
 PlotLabel -> "Fixed Charge Distribution (mol/m3)"]
Plot[cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1), {xStar, 0, 
  xStarMemMid}, PlotRange -> All, 
 PlotLabel -> "Coion Initial Condition (mol/m3)"]
Plot[-z2/z1*cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1), {xStar, 0,
   xStarMemMid}, PlotRange -> All, 
 PlotLabel -> "Counterion Initial Condition (mol/m3)"]

(*NDSolve*)
totalSol = 
 NDSolve[{c1Trnsprt, c2Trnsprt, elctrcTrnsprt, c1fIC, c1fDirichBC, 
   c2fIC, c2fDirichBC, phiDirichBC}, {c1f, c2f, phi}, {tStar, 0, 
   tRef}, {xStar, 0, xStarMemMid}(*{xStar}\[Element]xStarMesh*)(*,
  Method->{"PDEDiscretization"->{"FiniteElement",
  "MeshOptions"->{"IncludePoints"->xStarMesh["Coordinates"]}}}*), 
  Method -> {"PDEDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.002}}}]

With the solving method shown above, the solution does not seem correct. Therefore, I want to implement self-prescribed mesh points (similar to this post) to tackle region with sharp numerical change. My code to define the mesh point is:

(*Mesh Points*)
mshSze1 = 0.05;
mshSze2 = 0.005;
m = 15; (*Steepness*)
d[xStar_] = (mshSze1 - mshSze2 )/
     2*(Tanh[-m*(xStar - xStarMem0)] + 1) + mshSze2;(*Without dip*)
Plot[d[x], {x, 0, xStarMemMid}, PlotRange -> All, 
 PlotLabel -> "Mesh Size VS xStar"]

xlist = {0};
xend = xlist[[Length[xlist]]];
coord = {{1}};
i = 1;
While[xend < xStarMemMid,
 i = i + 1;
 coord = Join[coord, {{i}}];
 xend = xend + d[xend];
 If[xend > xStarMemMid, xend = xStarMemMid;, xend = xend;];
 xlist = Insert[xlist, xend, (Length[xlist] + 1)];
 ]
Needs["NDSolve`FEM`"]
xStarMesh = ToElementMesh[Map[{#} &, xlist]];
MeshRegion[xStarMesh]

And the NDSolve line changes to:

totalSol = 
 NDSolve[{c1Trnsprt, c2Trnsprt, elctrcTrnsprt, c1fIC, c1fDirichBC, 
   c2fIC, c2fDirichBC, phiDirichBC}, {c1f, c2f, phi}, {tStar, 0, 
   tRef}, {xStar}\[Element] xStarMesh] 

When I run the line above, the following errors appear:

NDSolve::femnotime: The differential equation cannot be solved as a time dependent equation as specified, most likely because the initial conditions given at (tStar==0.`) are not sufficient to define an initial value problem. As a consequence the differential equation will be solved as a time independent equation.

NDSolve::femimvr: The number of independent variables 2 ({xStar,tStar}) does not match the embedding dimension 1