My attempt to NDSolve a 2nd order nonlinear ODE
rmin = 10^(-40); (* as close to 0 as possible*)
bc = {u'[rmin] == 0, u'[1] == -u[1]};
ode = r*u''[r] + 2*u'[r] + r*(Pi/64)*Exp[-128*r^2]*(1 - u[r]^5) == 0;
s = NDSolve[{ode, bc}, u, {r, rmin, 1}, WorkingPrecision -> 70,
AccuracyGoal -> 20];
resulted in an accuracy of about $10^{-12}$ as can be seen from the plot
Plot[ode[[1]] /. s, {r, rmin, 1}, PlotRange -> All]
An alternative approach ( from this post)
ClearAll[s, u, v, rmin]
rmin = 10^(-40);
defv = u'[r] == v[r];
odev = r*v'[r] + 2 v[r] + r*(Pi/64)*Exp[-128*r^2]*(1 - u[r]^5) == 0;
bcv = {v[rmin] == 0, u[1] == -v[1]};
AbsoluteTiming[
s = NDSolve[{defv, odev, bcv}, {u, v}, {r, rmin, 1},
StartingStepSize -> 1*^-8, MaxStepSize -> 1*^-4,
PrecisionGoal -> 33, AccuracyGoal -> 33, WorkingPrecision -> 70,
MaxSteps -> 2*^5, InterpolationOrder -> All];]
resulted in an increased accuracy $10^{-19}$
Plot[odev[[1]] /. s, {r, rmin, 1}, PlotRange -> All
However computation time also increased significantly (about 1 min).
Is this the best accuracy Mathematica can achieve?
Motivation: The above boundary value problem corresponds to the construction of initial data for General Relativity free evolution (see Okawa, Cardoso, Pani, Phys.Rev.D, 90, 104032 (2014), eq.23 ).
Here is an attempt to address the full problem.
This kind of evolution ivolves only PDE's with respect to time and can in princple be performed in Mathematica.
However Mathematica probably lacks some tips and tricks used by Numerical GR algorithms in which case it will be more prone to crash due to numerical error on the initial data surface.
Method
used for integrating your ODE; I've found that for really stringent precision requirements,Method -> 'StiffnessSwitching"
(which uses Bulirsch-Stoer) does remarkably well. $\endgroup$