How to Solve PDE using Inactive Form

I am trying to solve this PDE using the inactive form, however I ran into error stating that the: The PDE coefficient 0. does not evaluate to a numeric matrix of dimensions {2,1} at the coordinate {0.5,0.5}; it evaluated to 0. instead. This error appears after using the new diffusion expression, diff2s which is essential for my code. Here is the code

ClearAll["Global`*"]

(*Parameters*)
Dcoeff = 9.24*10^-11; Xmax = 0.015; tmax = 10^8; A = 
 Dcoeff*tmax/Xmax^2;
Na = 6.022*10^23; vw = 2.99*10^-29; B1 = 50; D1 = 0.1;
CSinitial = 1; CFinitial = 1/1000; DelU = -5; vfsb = 
 6.5*10^-26; a = 2.5;

(*Binding site profile*)
cs[xstar_] := CSinitial*(0.5 - 0.5*Tanh[B1*(xstar - D1)]);

(*Bound cargo*)
cb[xstar_, CF_] := (Na*cs[xstar]*CF*vw)/(Exp[DelU] + Na*CF*vw);

(*Steric potential*)
Uster[CF_, xstar_] := -(1 + a)^2*
   Log[1 - cb[xstar, CF]*vfsb - CF*vfsb - cs[xstar]*vfsb];

(*Binding potential*)
Ubind := DelU;

(*Effective potential*)
Ueff[CF_, xstar_] := Uster[CF, xstar] + Ubind;

(*Diffusivity spatial function*)
phi[CF_, xstar_] := 1 - vfsb*(cb[xstar, CF] + CF + cs[xstar]);
Diff2s[CF_, xstar_] := (phi[CF, xstar]/(2 - phi[CF, xstar]))^2;


(*PDE in Inactive form*)
cEqn := (1 + (CSinitial/CFinitial)*
       D[cb[xstar, cfree[tstar, xstar]], cfree[tstar, xstar]])*
    D[cfree[tstar, xstar], tstar] == 
   A*Inactive[Div][{0, 
       Diff2s[cfree[tstar, xstar], 
         xstar]*(cfree[tstar, xstar]*
           Inactive[Grad][
             Ueff[cfree[tstar, xstar], xstar], {tstar, xstar}][[2]] + 
          Inactive[Grad][
            cfree[tstar, xstar], {tstar, xstar}][[2]])}, {tstar, 
       xstar}] + NeumannValue[0, xstar == 0];

(*Initial& boundary conditions*)
ic = cfree[0, xstar] == (0.5 - 0.5*Tanh[B1*(xstar - D1)]);
bc = DirichletCondition[cfree[tstar, xstar] == 0, xstar == 1];

(*Solve*)
sol = NDSolveValue[{cEqn, ic, bc}, 
   cfree, {tstar, 0, 1}, {xstar, 0, 1}, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.00005}}}];