As one of the possibly solutions we can decrease accuracy and use explicit method as follows
ClearAll["Global`*"]
(*Parameters*)
xRef = 10^(-8);
rhoRef = 0.1*nA;
j = 0.5; (*normalized conc/(m^2s)*)
rhoFStarVal = 120;(*normalized fixed charge concentration*)
rhoBStarVal = 20;(*normalized bulk concentration*)
phiBC = -0.25;(*in thermal voltage unit*)
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);
xL = 0;(*in xRef*)
xMemL = 0.75; (*in xRef*)
xMemR = 1.25;(*in xRef*)
xR = 2;(*in xRef*)
stepFunc[x_] = 0.5*(Tanh[100*x] + 1);
filterFunc[x_] = stepFunc[x - xMemL]*stepFunc[-(x - xMemR)];
rhoFStar[x_] = rhoFStarVal*filterFunc[x];
zf = -1;
z1 = 1;(*Counter ion valence*)
z2 = -1;(*Co ion valence*)
deltaU = -1/beta;
vW = 30*10^(-30);
debye2 = (kB*absTemp*eps*eps0)/(e^2*rhoRef);
(*Equations*)
elctrcTrnsprtV2 =
debye2/xRef^2*
phi''[x] == -(zf*rhoFStar[x] + z1*rho1fStar[x] + z2*rho2fStar[x])
ctionTrnsprt =
rho1fStar'[x] == -z1*rho1fStar[x]*phi'[x] - (j*xRef)/dfsvt
coionTrnsprt = rho2fStar'[x] == -z2*rho2fStar[x]*phi'[x]
(*NDSolve*)
totalSolV2 =
NDSolve[{ctionTrnsprt, coionTrnsprt, elctrcTrnsprtV2, phi[xL] == 0,
phi[xR] == phiBC, rho1fStar[xL] == rhoBStarVal,
rho2fStar[xL] == rhoBStarVal}, {phi, rho1fStar, rho2fStar}, {x, xL,
xR}, AccuracyGoal -> 5, PrecisionGoal -> 5,
Method -> "ExplicitEuler"]
(*Define the concentration function concentration[r,t]*)
phiSol = totalSolV2[[1]][[1]][[2]];
rho1fStarSol = totalSolV2[[1]][[2]][[2]];
rho2fStarSol = totalSolV2[[1]][[3]][[2]];
(*Solution Plotting*)
{Plot[phiSol[x], {x, xL, xR}, PlotRange -> All,
PlotLabel -> "phi (KbT/e), steady"],
Plot[rho1fStarSol[x], {x, xL, xR}, PlotRange -> {0, All},
PlotLabel -> "Free Counter Ion (0.1 mM/L), steady"],
Plot[rho2fStarSol[x], {x, xL, xR}, PlotRange -> Full,
PlotLabel -> "Free Co Ion (0.1 mM/L), steady"]}
