How to avoid NDSolve::ndsz problem (singularity problem)
Update notice 2: Found starting initial conditions that work in V10.
Update notice: bbgodfrey pointed out that the built-in shooting method does not work in V10.0.1 (nor V10.0.2).
Set up the shooting method by hand, so you can control the convergence. The problem is that small changes in the initial conditions cause the solutions to blow up before r reaches R. This confuses FindRoot, which is being used to solve for the initial conditions.
There are two ways to proceed in V9; in V10 only the first method works. (1) Do it completely by hand with ParametricNDSolve or (2) use the "Shooting" method option "StartingInitialConditions". We get singularity warnings with both approaches when FindRoot happens to try the wrong initial conditions. One can use Quiet to suppress these messages, or live with them. I used Quiet below, but normally I just live with them.
ParametricNDSolve (V10/V9)
sol = ParametricNDSolve[{eqs[[All, 1 ;; 2]], x1[R0] == a,
x2[R0] == b}, {x1, x2}, {r, R0, R}, {a, b},
"ExtrapolationHandler" -> {Automatic, "WarningMessage" -> False}];
bcs = Quiet[
FindRoot[{x1[a, b][R] == Cb, x2[a, b][R] == Sp} /. sol, {a,
0.1}, {b, 12}],
ParametricNDSolve::ndsz
]
(* {a -> 0.0000188874, b -> 20.8572} *)
{x1i, x2i} = {x1[a, b], x2[a, b]} /. bcs /. sol;
{Cb, Sp}
{x1i[R], x2i[R]}
(*
{8, 30}
{8., 30.}
*)
"StartingInitialConditions" (V9)
sol = Quiet[
NDSolve[eqs, {x1, x2}, {r, R0, R},
Method -> {"Shooting",
"StartingInitialConditions" -> {x1[R0] == 0.1, x1'[R0] == 0,
x2[R0] == 12, x2'[R0] == 0}}],
NDSolve::ndsz
];
{Cb, Sp}
{x1[R], x2[R]} /. sol
(*
{8, 30}
{{8.00087, 30.001}}
*)
Update: "StartingInitialConditions" (V10)
The problem is that the solution {a -> 0.0000188807, b -> 20.8572} is close to the edge of the domain where solution is nonsingular for r < R = 0.4. The use of FindRoot in the built-in shooting method tends to step outside this domain if the starting value for x1[R0] == a is to large. In V10 (in this example at least), this seems to cause the shooting method to fail. The solution is to start with a value of a that is between the solution and the boundary of the domain. Instead of figuring all that out, it seems easier to use the first approach above, though.
sol = NDSolve[eqs, {x1, x2}, {r, R0, R},
Method -> {"Shooting",
"StartingInitialConditions" -> {x1[R0] == 0, x1'[R0] == 0,
x2[R0] == 20., x2'[R0] == 0}}]
The solution is (virtually) the same as the V9 solution.
I think part of the issue lies with the coordinate singularity at r=0. I made the following transformations, first x1 -> x1/KO2 , x2->x2/KS . Then u1=x1/r and u2=x2/r. I also assumed that you were solving from R0>0 to R because you couldn't solve at R0=0. I set the boundary conditions u1[0]=0, u2[0]=0 . These transformations were used because
Laplacian[ f[r]/r, {r, θ, z}, "Spherical"] // Expand
(* = f''[r]/r *)
I found solutions for R=0.7
Then I wrote the equations as:
Clear[R, px, Mm, Ds, DO2, YXS, KS, Sp, Cb, KO2, R0, YXO, r, u1, u2, eqs2];
A = (Mm px)/(YXS DO2 KO2);
B = (Mm px)/(YXO Ds KS);
V = (r^3 u1[r])/(1 + r u1[r]) u2[r]/(1 + r u2[r]);
Cbhat = Cb/KO2;
Sphat = Sp/KS;
eqs2 = {u1''[r] == A V, u1[0] == 0, u1[R] == R Cbhat, u2''[r] == B V,
u2[0] == 0, u2[R] == R Sphat}
R = .7; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.1; R0 = 0.000; YXO = 0.42857;
sol = NDSolve[eqs2, {u1, u2}, {r, 0, R}]
Plot[{(Evaluate[u1[r]] /. sol), R Cbhat}, {r, R0, R}, PlotRange -> {0, R Cbhat}]
Plot[{(Evaluate[u2[r]] /. sol), R Sphat}, {r, R0, R}, PlotRange -> {0, R Sphat}]
Beyond this I would try the shooting method suggested above. Beyond that I would try to scale u1, u2, and r to get rid of the very large parameter A. With the proper scaling a perturbative approach will suggest itself. If the solution actually ceases to exist for some value of R then an analytic approach will help. For starters, solutions disappear in pairs and at the bifurcation point the two solutions are equal. Thus, if you can find two solutions that look as if they're becoming equal as you vary R it would be indicative of a "saddle-node" bifurcation (of the sort that occurs at a=0 for the equation a-x^2=0).
Afterthoughts:
Clearly there is a boundary layer occurring for larger R. The fact that a numerical method fails to return a solution there is not surprising, nor is it necessarily indicative of a bifurcation point. There are various other tricks you might consider. This looks to be some sort of chemistry problem for which you're looking for stationary solutions? Either way you might consider solving the time-dependent version of the problem. It will always have a solution and if your steady solution is stable then the time dependent version will likely settle to it, assuming you have initial conditions that are in the basin of attraction for that solution.
COMMENT: MAJOR ERRORS CORRECTED IN THE BELOW and analytic solution broached
I should mention that you can reduce the dimension of the boundary value problem by noting that B u1"-A u2" = 0, so that w=B u1-A u2 obeys w" = 0 -> w = alpha r +beta and the boundary condtions give you beta=0 and alpha = B Cbhat - A Sphat. So you can keep the equation for u1 and replace u2= B/A u1 - (B Cbhat - A Sphat)/A . If you do this then you'll get the following equations:
Clear[R, px, Mm, Ds, DO2, YXS, KS, Sp, Cb, KO2, R0, YXO, r, u1, u2, eqs2];
A = (Mm px)/(YXS DO2 KO2)
B = (Mm px)/(YXO Ds KS);
V = (r^3 u1[r])/(1 + r u1[r]) u2[r]/(1 + r u2[r]);
V = V /. u2[r] -> ( B/A u1[r] - (B Cbhat - A Sphat) r/A)
(*
V = V/.u1[r]\[Rule]((Cbhat+Sphat)r+A/B u2[r])
*)
Cbhat = Cb/KO2;
Sphat = Sp/KS;
eqs2 = {u1''[r] == A V, u1[0] == 0, u1[R] == R Cbhat}
R = .8; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.1; R0 = 0.000; YXO = 0.42857;
A = 10^8
B
sol = NDSolve[eqs2, u1, {r, 0, R}];
u2[r_] = B/A ( Evaluate[u1[r]] /. sol) - (B Cbhat - A Sphat) r/A;
Plot[{(Evaluate[u1[r]] /. sol), R Cbhat}, {r, R0, R}, PlotRange -> {0, R Cbhat}, PlotStyle -> {Black, Red, Blue}]
Plot[{u2[r], Sphat r}, {r, 0, R}, PlotStyle -> {Black, Red}]
Note that I've overridden the defined value of A and set A=10^8. It wouldn't have mattered if I'd set it to 10^google. This is asymptopia for A. If you assume that u1 is bounded as A-> infinity (it is) then for u2 you get u2 = Sphat r in the large A limit. This still leaves you needing to solve:
u1" = A V but now, u2-> Sphat r in V. You'll have to figure out what the u2 solution means means but it means something in terms of the physics of your problem. If you pursue this you'd probably want to compute a correction for u2 as A=14,000 is not yet asymptotic. Now there is clearly a boundary layer and it is clearer but there's still work to be done to tease out the physics.
Below is the code that solves the 1-variable problem that remains after setting u2=Sphat r. I'll post this now but I think there are a few more easy steps to make.
Clear[R, px, Mm, Ds, DO2, YXS, KS, Sp, Cb, KO2, R0, YXO, r, u1, u2, eqs2];
A = (Mm px)/(YXS DO2 KO2)
B = (Mm px)/(YXO Ds KS);
V = (r^3 u1[r])/(1 + r u1[r]) u2[r]/(1 + r u2[r]);
V = V /. u2[r] -> ( Sphat r)
Cbhat = Cb/KO2;
Sphat = Sp/KS;
eqs2 = {u1''[r] == A V, u1[0] == 0, u1[R] == R Cbhat}
R = 1; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.1; R0 = 0.000; YXO = 0.42857;
A = 100000000
B
sol = NDSolve[eqs2, u1, {r, 0, R}];
u2[r_] = B/A ( Evaluate[u1[r]] /. sol) - (B Cbhat - A Sphat) r/A;
LogLogPlot[{(Evaluate[u1[r]] /. sol), R Cbhat}, {r, .00000001, R}, PlotRange -> {0, R Cbhat}, PlotStyle -> {Black, Red, Blue}]
Plot[{u2[r], Sphat r}, {r, 0, R}, PlotStyle -> {Black, Red}]
Anther approach in V12.0 is to use the finite element method. For that we remove the zero derivatives (these are zero NeumannValues)
eqs = {{x1''[r] + 2/r x1'[r] == Vo/DO2,
x1[R] == Cb}, {x2''[r] + 2/r x2'[r] == Vs/Ds, x2[R] == Sp}}
with the remaining parameters:
Vs = px 1/YXS (Mm x2[r])/(KS + x2[r]) x1[r]/(KO2 + x1[r]);
Vo = px 1/YXO (Mm x2[r])/(KS + x2[r]) x1[r]/(KO2 + x1[r]);
R = 0.4; px = 32000; Mm = 0.1; Ds = 9; DO2 = 7.2; YXS = 0.3; KS = 10;
Sp = 30; Cb = 8; KO2 = 0.2; R0 = 0.000001; YXO = 0.42857;
We call NDSolveValue
sol = NDSolveValue[eqs, {x1, x2}, {r, R0, R},
Method -> "FiniteElement"];
Looking at the solution:
Through[sol[R]]
{8., 30.}
This what the other solutions found.