Problem Encountered when Solving a System Consisting of Two PDEs and an ODE in a Semi-NDSolve-based Approach

This is a continuation for question 300522.

Firstly, I have changed the equation for $$C_{2\text{b}}^*$$ in my system. The updated $$C_{2\text{b}}^*$$ is expressed as below:

$$$$C_{2\text{b}}^* = C_{F}^* \frac{C_{2\text{f}} N_A v_w}{C_{2\text{f}} N_A v_w +\exp{\beta \Delta U} }$$$$

For a clearer context, I am basically solving for a transient Poisson-Nernst-Planck (PNP) system. It consists of three major equations: one transport PDE for c1f (co-ion), one transport PDE for c2f (counter-ion), and one poisson equation (the ODE) that governs the electrical potential phi induced by the salt ion and membrane fixed charge group. Additionally, I also consider the binding of the counter-ion at the membrane fixed charge, termed as c2b or $$C_{2\text{b}}^*$$.

The computational domain consists of the salt solution (0<xStar<3) and the membrane (3<xStar<6). Here, I consider one property difference between the salt solution and the membrane, which is the existence of fixed-charge cF[xStar] in the membrane, while none in the salt solution. This is a piecewise function that is approximated by xStar-dependent hyperbolic tangent function. To note, it is expected to have a steep change in mathematical solution value at the solution-membrane interface, and I believe this might be one of the possible origins for the solving issue.

I believe that a logically-sound numerical solution should exist for my system. This is because I managed to obtain the solution at tStar=0 and steady state.

My code to obtain solution at tStar=0:

(*Phi Solution for Initial Condition*)
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 = 6;
cSaltBar = 1500;(*mol/m3*)
cFVal = 300;(*mol/m3*)
cF[xStar_] = cFVal/2*(Tanh[15*(xStar - xStarMem0)] + 1);
z1 = -1;
z2 = 1;
deltaU = -1/beta;
vW = 30*10^(-30);

(*BCs*)
phiDirichBC = phi[0] == 0;
phiNuemBC = NeumannValue[0, xStar == xStarMemMid];

(*IC Equations*)
c1fIC[xStar_] = -z2/z1*cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1);
c2fIC[xStar_] = cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1);
c2bIC[xStar_] =
cF[xStar]*(nA*c2fIC[xStar]*vW)/(nA*c2fIC[xStar]*vW +
Exp[beta*deltaU]);
elctrcTrnsprtIC =
D[eps*eps0/(xRef^2*e*nA)*D[phi[xStar], xStar],
xStar] == -(-cF[xStar] + z2*(c2fIC[xStar] + c2bIC[xStar]) +
z1*c1fIC[xStar]) + phiNuemBC;

(*NDSolve*)
totalSol =
NDSolve[{elctrcTrnsprtIC, phiDirichBC}, {phi}, {xStar, 0,
xStarMemMid},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" ->
0.0001}}}];
(*Define the concentration function concentration[r,t]*)
phiSolIC[xStar_] = phi[xStar] /. totalSol[[1]];

(*Solution Plotting*)
Plot[Evaluate@phiSolIC[xStar], {xStar, 0, xStarMemMid},
PlotRange -> All, FrameLabel -> {"xStar", "Elec Pot (V)"},
PlotLabel -> "Electrical Potential, Initial"]
Plot[c1fIC[xStar], {xStar, 0, xStarMemMid}, PlotRange -> {0, Full},
FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Free Co Ion, Initial"]
Plot[c2fIC[xStar], {xStar, 0, xStarMemMid}, PlotRange -> {0, Full},
FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Free Ct Ion, Initial"]
Plot[c2bIC[xStar], {xStar, 0, xStarMemMid}, PlotRange -> {0, Full},
FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Binded Ct Ion, Initial"]


My code to obtain steady-state solution:

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 = 6;
cSaltBar = 1500;(*mol/m3*)
cFVal = 300;(*mol/m3*)
cF[xStar_] = cFVal/2*(Tanh[15*(xStar - xStarMem0)] + 1);
z1 = -1;
z2 = 1;
deltaU = -1/beta;
vW = 30*10^(-30);

(*BCs*)
phiDirichBC = phi[0] == 0;
phiNuemBC = NeumannValue[0, xStar == xStarMemMid];

(*Equations*)
c1f[xStar_] = -z2/z1*cSaltBar*Exp[-beta*z1*e*phi[xStar]] ;
c2f[xStar_] = cSaltBar*Exp[-beta*z2*e*phi[xStar]] ;
c2b[xStar_] =
cF[xStar]*(nA*c2f[xStar]*vW)/(nA*c2f[xStar]*vW + Exp[beta*deltaU]);
elctrcTrnsprt =
D[eps*eps0/(xRef^2*e*nA)*D[phi[xStar], xStar],
xStar] == -(-cF[xStar] + z2*(c2f[xStar] + c2b[xStar]) +
z1*c1f[xStar]) + phiNuemBC;

(*Input Parameter Plotting*)
Plot[cF[x], {x, 0, xStarMemMid}, PlotRange -> All,
PlotLabel -> "Fixed Charge Distribution (mol/m3)"]

(*NDSolve*)
totalSol =
NDSolve[{elctrcTrnsprt, phiDirichBC}, {phi}, {xStar, 0,
xStarMemMid},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}}];

(*Define the concentration function concentration[r,t]*)
phiSol[xStar_] = phi[xStar] /. totalSol[[1]];

(*Solution Plotting*)
Plot[Evaluate@phiSol[xStar], {xStar, 0, xStarMemMid},
PlotRange -> All, FrameLabel -> {"xStar", "Elec Pot (V)"},
Plot[c1f[xStar] /. {phi -> phiSol}, {xStar, 0, xStarMemMid},
PlotRange -> {0, All}, FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Free Co Ion, steady"]
Plot[c2f[xStar] /. {phi -> phiSol}, {xStar, 0, xStarMemMid},
PlotRange -> {0, All}, FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Free Ct Ion, steady"]
Plot[c2b[xStar] /. {phi -> phiSol}, {xStar, 0, xStarMemMid},
PlotRange -> {0, All}, FrameLabel -> {"xStar", "Conc (mol/m3)"},
PlotLabel -> "Binded Ct Ion, steady"]


To solve the Poisson-Nernst-Planck in transient state, I have adopted the Method of Lines approach from question 78493, which is commonly used to solve for integro-differential equation and PDE-ODE-coupled system. To discretize my system, I am using the pdetoode function from the accepted answer of question 127980. Not only that, I also differentiated both sides of my ODE wrt timeStar, making it a PDE. Subsequently, I included the initial condition for phi, which is the variable to be solved from my ODE.

Here is the code for pdetoode:

ClearAll["Global*"]
Clear[fdd, pdetoode, tooderule, pdetoae, diffbc, rebuild]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolveFiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]],
rest];
pdetoode[{func__}[var__], rest__] :=
pdetoode[Alternatives[func][var], rest];
pdetoode[front__, grid_?VectorQ, o_Integer, periodic_ : False] :=
pdetoode[front, {grid}, o, periodic];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer,
periodic : True | False | {(True | False) ..} : False] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
Piecewise[{{1, PossibleZeroQ[# - #2[[1]]]}, {-1,
PossibleZeroQ[# - #2[[-1]]]}}, All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};

fdd[{dx1, dx2}, {grid},
Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o,
PeriodicInterpolation -> periodic]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]

pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] :=
Module[{t},
Function[
pde, #[pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]]

diffbc[rst__][a : _List | _Equal] := diffbc[rst] /@ a
diffbc[dvar : {t_, order_} | (t_) .., sf_ : 0][a_] /; sf =!= t :=
sf a + D[a, dvar]

rebuild[funcarray_, grid_?VectorQ, timeposition_ : 1] :=
rebuild[funcarray, {grid}, timeposition]

rebuild[funcarray_, grid_, timeposition_?Negative] :=
rebuild[funcarray, grid, Range[Length@grid + 1][[timeposition]]]

rebuild[funcarray_, grid_, timeposition_ : 1] /;
Dimensions@funcarray === Length /@ grid :=
With[{depth = Length@grid},
ListInterpolation[
Transpose[
Map[DeveloperToPackedArray@#["ValuesOnGrid"] &, #, {depth}],
Append[Delete[Range[depth + 1], timeposition], timeposition]],
Insert[grid, Flatten[#][[1]]["Coordinates"][[1]],
timeposition]] &@funcarray]


Below is my code:

(*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 = 6;
cSaltBar = 1500;(*mol/m3*)
cFVal = 300;(*mol/m3*)
cF[xStar_] = cFVal/2*(Tanh[15*(xStar - xStarMem0)] + 1);
z1 = -1;
z2 = 1;
deltaU = -1/beta;
vW = 30*10^(-30);

(*Governing Equations*)
c2b[xStar_, tStar_] =
cF[xStar]*(nA*c2f[xStar, tStar]*vW)/(nA*c2f[xStar, tStar]*vW +
Exp[beta*deltaU]);
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];
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];
elctrcTrnsprt =
D[D[eps*eps0/(xRef^2*e*nA)*D[phi[xStar, tStar], xStar], xStar],
tStar] == -D[(-cF[xStar] +
z2*(c2f[xStar, tStar] + c2b[xStar, tStar]) +
z1*c1f[xStar, tStar]), tStar];

(*Phi Solution for Initial Condition*)
(*IC Equations*)
c1fICsol[xStar_] = -z2/z1*
cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1);
c2fICsol[xStar_] = cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1);
c2bICsol[xStar_] =
cF[xStar]*(nA*c2fICsol[xStar]*vW)/(nA*c2fICsol[xStar]*vW +
Exp[beta*deltaU]);
elctrcTrnsprtIC =
D[eps*eps0/(xRef^2*e*nA)*D[phiIC[xStar], xStar],
xStar] == -(-cF[xStar] + z2*(c2fICsol[xStar] + c2bICsol[xStar]) +
z1*c1fICsol[xStar]) + phiICNuemBC;
(*BCs*)
phiICDirichBC = phiIC[0] == 0;
phiICNuemBC = NeumannValue[0, xStar == xStarMemMid];
(*NDSolve*)
totalSol =
NDSolve[{elctrcTrnsprtIC, phiICDirichBC}, {phiIC}, {xStar, 0,
xStarMemMid},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" ->
0.0001}}}(*,{xStar}\[Element]xStarMesh*)];
(*Define the concentration function concentration[r,t]*)
phiSolIC[xStar_] = phiIC[xStar] /. totalSol[[1]]

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

(*Spatial Discretization of the Governing Eqns, BCs, and ICs*)
points = 25; (*Intervals is points - 1*)
domain = {0, xStarMemMid};
grid = Array[# &, points, domain];
difforder = 1;
ptoofunc =
pdetoode[{c1f, c2f, phi}[xStar, tStar], tStar, grid, difforder];
removeredundant = #[[2 ;; -2]] &;
dscrtzdBCs = {c1fDirichBC, c1fNuemBC, c2fDirichBC, c2fNuemBC,
phiDirichBC, phiNuemBC} // ptoofunc;(*//MatrixForm*)
dscrtzdICs =
Map[removeredundant,
ptoofunc[{c1fIC, c2fIC, phiIC}]];(*///MatrixForm*)
dscrtzdGovEqns =
Map[removeredundant,
ptoofunc[{c1Trnsprt, c2Trnsprt, elctrcTrnsprt}]];(*//MatrixForm*)

(*Solving*)
sollst =
NDSolveValue[{dscrtzdBCs, dscrtzdICs, dscrtzdGovEqns},
Outer[#[#2] &, {c1f, c2f, phi}, grid], {tStar, 0, 0.1},
MaxSteps -> Infinity];
(*Rebuild the solution for the PDE from the solution for the ODE set:*)
sol = rebuild[sollst, grid];


However, I encountered the following issue despite I have equal number of discretized governing equations and discretized initial conditions.

NDSolveValue::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

Code to check my number of discretized governing equations/discretized initial conditions/discretized boundary conditions:

Length@Flatten@dscrtzdGovEqns
Length@Flatten@dscrtzdICs
Length@Flatten@dscrtzdBCs


I am also curious to know if my system can be solved entirely just by using NDSolve, without the semi-manually-done spatial discretization.

• Currently you're discretizing the system to a DAE system, which is harder to solve compared with ODE system. (Notice the diffbc function in the package, search in this site and you'll see example using it. ) But even if I improve the definition of dscrtzdBCs and dscrtzdICs, NDSolve is having difficulty with the system. According to my limited experience, this can be a sign of wrong equation or improper time interval. So, let me ask: 1. Are you sure the current end time i.e. 0.1 is proper? 2. Are you sure the equation system is correct? Commented Mar 22 at 8:00
• @xzczd Actually, I am not sure of the end time. I am still estimating it now, thats why I simulate it from 0 to 0.01 seconds (which is equivalent to 0.1 in tStar domain, which is the scaled time). For a clearer context, I am computing the transport of salt ions (one PDE for c1f which is the co-ion, and one PDE for c2f which is the counter-ion) in a domain which consists of the solution (0<xStar<3) and the membrane (3<xStar<6). For completeness, I also consider poisson equation (the ODE) that governs the electrical charges caused by the salt ion and fixed charge group in the membrane. Commented Mar 22 at 14:36
• Oh the domain is consist of different materials? Then that'll indeed be an issue. You need to avoid the automatic differentiation of the outer D, see e.g. the discussion here. Commented Mar 22 at 14:37
• One important thing to note is that, there is a steep change in mathematical solutions between the membrane and the contacting salt solution. I suspect this is one of the origins of the issue too. Commented Mar 22 at 14:37
• @xzczd the domain consists of two different regimes (with different materials, so different transport properties). To tackle this, I am using a hyperbolic tangent function to simulate the jump in values of properties. Currently, the only difference in the properties that i consider is the fixed charge, cF[xStar]. It exists in the membrane (3<xStar<6), but disappears in the contacting salt solution (0<xStar<3). This transition is approximated by a hyperbolic tangent function. Commented Mar 22 at 14:41

This problem can be solved with using implicit Euler for time step and linear FEM for space. This method described, for example, here and here. First we transform equations to the system of ODEs and solve it with NDSolve as follows

(*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 = 6;
cSaltBar = 1500;(*mol/m3*)
cFVal = 300;(*mol/m3*)
cF[xStar_] = cFVal/2*(Tanh[15*(xStar - xStarMem0)] + 1);
z1 = -1;
z2 = 1;
deltaU = -1/beta;
vW = 30*10^(-30);
mu = 10^-2; nu = eps*eps0/(xRef^2*e*nA);
(*Governing Equations*)
c2b[xStar_, tStar_] =
cF[xStar]*(nA*c2f[xStar, tStar]*vW)/(nA*c2f[xStar, tStar]*vW +
Exp[beta*deltaU]);
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];
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];
elctrcTrnsprt =
mu D[phi[xStar, tStar], tStar] ==
D[D[phi[xStar, tStar], xStar],
xStar] + (-cF[xStar] +
z2*(c2f[xStar, tStar] + c2b[xStar, tStar]) +
z1*c1f[xStar, tStar])/nu;

(*Phi Solution for Initial Condition*)
(*IC Equations*)
c1fICsol[xStar_] = -z2/z1*
cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1);
c2fICsol[xStar_] = cSaltBar/2*(Tanh[-15*(xStar - xStarMem0)] + 1);
c2bICsol[xStar_] =
cF[xStar]*(nA*c2fICsol[xStar]*vW)/(nA*c2fICsol[xStar]*vW +
Exp[beta*deltaU]);
elctrcTrnsprtIC =
D[eps*eps0/(xRef^2*e*nA)*D[phiIC[xStar], xStar],
xStar] == -(-cF[xStar] + z2*(c2fICsol[xStar] + c2bICsol[xStar]) +
z1*c1fICsol[xStar]) + phiICNuemBC;
(*BCs*)
phiICDirichBC = phiIC[0] == 0;
phiICNuemBC = NeumannValue[0, xStar == xStarMemMid];
(*NDSolve*)
{phiSolIC} =
NDSolveValue[{elctrcTrnsprtIC, phiICDirichBC}, {phiIC}, {xStar, 0,
xStarMemMid},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" ->
0.002}}}(*,{xStar}\[Element]xStarMesh*)];
(*Define the concentration function concentration[r,t]
phiSolIC[xStar_]=phiIC[xStar]/. totalSol[[1]];*)

mesh = phiSolIC["ElementMesh"]

Out[]= ElementMesh[{{0., 6.}}, {LineElement["<" 3000 ">"]}]

C1F[0][x_] := -z2/z1*cSaltBar/2*(Tanh[-15*(x - xStarMem0)] + 1);
C2F[0][x_] := cSaltBar/2*(Tanh[-15*(x - xStarMem0)] + 1);
Phi[0][x_] := phiSolIC[x]; dt = 1/500;

AbsoluteTiming[Do[C2B = cF[x]*((nA*C2F[i - 1][x]*vW)/(nA*C2F[i - 1][x]*vW + Exp[beta*deltaU]));
{C1F[i], C2F[i], Phi[i]} = NDSolveValue[{(c1f[x] - C1F[i - 1][x])/dt ==
(dfsvt*tRef*(beta*e*z1*Derivative[1][c1f][x]*Derivative[1][Phi[i - 1]][x] + Derivative[2][c1f][x] +
beta*e*z1*C1F[i - 1][x]*Derivative[2][phi][x]))/xRef^2, (cF[x]*((nA*c2f[x]*vW)/(nA*C2F[i - 1][x]*vW + Exp[beta*deltaU])) + c2f[x] -
C2B - C2F[i - 1][x])/dt == (dfsvt*tRef*(beta*e*z2*Derivative[1][c2f][x]*Derivative[1][Phi[i - 1]][x] + Derivative[2][c2f][x] +
beta*e*z2*C2F[i - 1][x]*Derivative[2][phi][x]))/xRef^2, (eps*eps0*Derivative[2][phi][x])/(e*nA*xRef^2) ==
(-z1)*c1f[x] - z2*(cF[x]*((nA*c2f[x]*vW)/(nA*C2F[i - 1][x]*vW + Exp[beta*deltaU])) + c2f[x]) + cF[x],
DirichletCondition[{c1f[x] == (-z2/z1)*cSaltBar, c2f[x] == cSaltBar, phi[x] == 0}, x == 0]}, {c1f, c2f, phi}, Element[{x}, mesh]]; ,
{i, 1, 50}]]


It takes about 38s on my laptop. Visualization

{Plot[Evaluate[Table[C1F[i][x], {i, 0, 50, 5}]], Element[{x}, mesh],
PlotLegends -> Table[Row[{"t =", N[i  dt]}], {i, 0, 50, 5}],
PlotRange -> All, Frame -> True, FrameLabel -> {"x", "c1f"}],
Plot[Evaluate[Table[C2F[i][x], {i, 0, 50, 5}]], Element[{x}, mesh],
PlotLegends -> Table[Row[{"t =", N[i  dt]}], {i, 0, 50, 5}],
PlotRange -> All, Frame -> True, FrameLabel -> {"x", "c2f"}]}



Visualization at point x=3.3

{ListLinePlot[Table[{dt i, C1F[i][3.3]}, {i, 0, 50}], Frame -> True,
FrameLabel -> {"t", "c1f"}, PlotRange -> All],
ListLinePlot[Table[{dt i, C2F[i][3.3]}, {i, 0, 50}], Frame -> True,
FrameLabel -> {"x", "c2f"}, PlotRange -> All],
ListLinePlot[Table[{dt i, Phi[i][3.3]}, {i, 1, 50}], Frame -> True,
FrameLabel -> {"t", "phi"}, PlotRange -> All]}


Electrical potential at different times (abs value)

LogPlot[Evaluate[Table[Abs[Phi[i][x]], {i, 0, 50, 5}]], {x, 0, 6},
PlotRange -> Automatic, Frame -> True, FrameLabel -> {"x", "Phi"},
PlotLegends -> Table[Row[{"t =", N[i  dt]}], {i, 0, 50, 5}]]


• Thank you so much for the idea! I just checked the solution. However, it does not seem to match the trend of the steady-state solution that I provided in the post. I made a change to the mesh points. My code to create mesh is shown in the following: Commented Mar 23 at 17:45
• (*Mesh Points*) mshSze1 = 0.05; mshSze2 = 0.005; m = 5; (*Steepness*) (*d[xStar_]=(mshSze1-mshSze2 \ )/2*(Tanh[-m*(xStar-xStarMem0)]+1)+mshSze2;(*Without dip*)*) d[xStar_] = (mshSze1 - mshSze2 )/ 2*(Tanh[-m*(xStar - (xStarMem0 - 0.5))] + 1) + (mshSze1 - mshSze2 )/ 2*(Tanh[m*(xStar - (xStarMem0 + 0.5))] + 1) + mshSze2; Plot[d[x], {x, 0, xStarMemMid}, PlotRange -> {mshSze2, mshSze1}, PlotLabel -> "Mesh Size VS xStar"] Commented Mar 23 at 17:45
• 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["NDSolveFEM"] mesh = ToElementMesh[Map[{#} &, xlist]]; Commented Mar 23 at 17:46
• With this change being made, I realized that the solution reaches steady state within the first few frames of the computation. As I further reduce dt ( sub-microsecond range, note that the actual time is i*dt*tRef, while the scaled time is i*dt), the solution begins to exhibit noticeable instabilities, before I am able to observe the detailed transient behavior. Commented Mar 23 at 17:53
• @Johnson No it is not a typo. This is simplification to force NDSolve to activate linear FEM, while with your correction NDSolve` uses nonlinear FEM. Commented Mar 24 at 22:17