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"]} [![Figure 1][1]][1] [1]: https://i.sstatic.net/kmSFkdb8.png