NDSolve ODE with Hyperbolic Trig Func: Neumann Boundary Values Unsatisfied

Question

I'm trying to solve a tricky 2nd order ODE (the 1D Poisson-Boltzmann equation) with NDSolve for the electric potential at a distance (r) from the center of a charged particle of radius (a) at a density (d) with charge (q); it has two Neumann boundary conditions (BCs). It's an intricate equation and it isn't surprising it's hard [for me] to get a numerical solution. I feel I am close, but as yet no cigar...

soln =
 With[{a = 56/2*10^-9, d = 2.22814*10^19, q = 900},
   NDSolve[
{f''[r] + (2 f'[r])/r == (1.03975*10^8)^2 Sinh[f[r]],
 f'[a] == -((q (.714295*10^-9))/a^2),
 f'[a*((4. π d a^3)/3.)^(-1/3)] == 0}, f,
{r, 2*10^-8, 2.5*10^-7}]] // Flatten

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations.>>

NDSolve::berr: The scaled boundary value residual error of 8.240798196770146`*^9 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found.>>

I feel the two are likely related - if it had more iterations it may converge to a solution that satisfies the boundary values. This explains the first of the following attempts at resolving the issue: (listed in the order I suspect might be most valuable)

• Instructing it to run more iterations, but MaxSteps doesn't do the job and MaxIterations is not supported as an option of NDSolve; I feel inept not finding how to specify this.

• I've read that using Method -> {"Shooting","StartingInitialConditions" -> {f'[a] ==1.25((q (.714295*10^-9))/a^2)}} i.e. simply adding a coefficient to start slightly (and I've tried not so slightly) away from the initial condition might work. This returns singularity / stiffness errors. This question seems to indicate that if I choose the starting IC more wisely, it may eliminate the stiffness issues, but I'm at a loss as to what else to try.

• Changing the coefficient of the RHS of the equation (the value 1.03975*10^8) corresponding to a parameter that will eventually be fed into ParametricNDSolve. This leads to stiffness issues that can be resolved with Method -> "StiffnessSwitching, only to return the same errors quoted above.

• Distances (a), (r), that in the density (d), etc. are in meters; I've similar problems in the more natural length scale of nanometers.

• Lowering the AccuracyGoal and PrecisionGoal. Simply increases the discrepancy with the BC.

• Adding understood BC at infinity (f[inf] = 0). Mathematica crashes, understandably.

Does anyone have any idea how to resolve the issue? For convenience, the function to plot the result of the integration over the concerning range (with a gridlines at each BC point, where the green ones are the value of the function f[r] itself that I really aim to acquire):

Module[{a = 56/2*10^-9, d = 2.22814*10^19},
 Plot[{f[r]} /. soln, {r, 2*10^-8, 2.5*10^-7}, PlotStyle -> Thick,
  GridLines -> {{{a, Red}, {2.20457*10^-7, Darker@Green}},  
   {{f[a*((4 π d a^3)/3.)^(-1/3)] /. soln, Darker@Green}}}]]

Answer 1

A step backward i.e. turn to finite difference method (FDM) seems not to be a bad choice here. I'll use pdetoae for the generation of difference equations:

With[{a = 56/(2 10^9), d = 2.22814 10^19, q = 900}, 
 domain = Rationalize[#, 0] &@{a, a ((4. π d a^3)/3.)^(-1/3)};
 eq = f''[r] + (2 Derivative[1][f][r])/r == (1.03975 10^8)^2 Sinh[f[r]];
 bc = Rationalize[#, 0] &@{f'[a] == -((q .714295)/(10^9 a^2)), 
    f'[a ((4. π d a^3)/3.)^(-1/3)] == 0};]

points = 100;
grid = Array[# &, points, domain];
difforder = 2;
(* Definition of pdetoae isn't included in this code piece,
   please find it in the link above. *)
ptoafunc = pdetoae[f[r], grid, difforder];
ae = ptoafunc[eq][[2 ;; -2]];
aebc = ptoafunc[bc];
soldata = FindRoot[Rationalize[{ae, aebc} // Flatten, 0], {f@#, 1/2} & /@ grid, 
    WorkingPrecision -> 16][[All, -1]];
ListLinePlot[soldata, DataRange -> domain, PlotRange -> All]
(*sol=ListInterpolation[soldata,domain];
Plot[sol[x],{x,domain[[1]],domain[[2]]}]*)