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}]
FindRooterrors -- ??? $\endgroup$