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

Question

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:

$$ \begin{equation} 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} } \end{equation} $$

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)"}, 
 PlotLabel -> "Electrical Potential, steady"]
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__] := NDSolve`FiniteDifferenceDerivative@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, 
        MapThread[
         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[Developer`ToPackedArray@#["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.

Answer 1

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}]]