I have a relatively simple 2nd Order ODE problem I'm trying to solve numerically

$p'' + \frac{1}{r}p' = A(\frac{p(r)}{p(r) + k_{m}})$

where $A = 7.5795*10^9$ and $k_{m} = 1$.

So trying to get numerical solution in Mathematica between the boundaries $r_{o} = 5 \times 10^{-6}$ and $r_{c}$ ($r_{c}$ is at LEAST 0.000102329 and likely 10-30 microns bigger, but I'll use small estimate for this example . We also know $p(r_{o}) = p_{o} = 100$, and $p'(r_{c}) = p(r_{c}) = 0$; but this is seriously upsetting Mathematica; I'm getting the following warnings depending on what versions I run;

s = NDSolve[{D[r*O2'[r], r] - con*r*((O2[r])/(O2[r] + km)) == 0,  O2[rc] == 0, O2[ro] == po}, O2, {r, ro, rc},Method -> {"StiffnessSwitching"}, MaxSteps -> 100000]

Produces a host of errors of the form;

NDSolve::ndsz: At r == 0.00009971646984096604`, step size is effectively zero; >singularity or stiff system suspected

Trying the other condition

s = NDSolve[{D[r*O2'[r], r] - con*r*((O2[r])/(O2[r] + km)) == 0,  O2'[rc] == 0, O2[ro] == po}, O2, {r, ro, rc},Method -> {"StiffnessSwitching"}, MaxSteps -> 100000]

produces the errors

FindRoot::cvmit: Failed to converge to the requested accuracy or precision within 100 >iterations. NDSolve::berr: "There are significant errors {-0.125305,-99.99} in the boundary value >residuals. Returning the best solution found."

The attempted solutions look frankly obscenely wrong, so Im guessing there's a stiffness or singularity probem. I've tried a few tricks to try and circumvent this and I'm getting no joy; is there something I'm missing? So far I've tried increasing the value for rc, and re-writing the equation differently but I'm still having problems. Anyone any idea of why I'm getting these errors and how to resolve them?

  • $\begingroup$ In my understanding, you'd expect a singularity at r = 0, but the function is not evaluated here, and there should be no singular points between ro and rc. Most odd... $\endgroup$
    – DRG
    Oct 21, 2013 at 14:59

1 Answer 1


Your equation can be solved by "shooting method" with a slight adjustment of rc:

con = 7.5795 10^9; km = 1; po = 100;
ro = 5 10^-6; rc = 132.329 10^-6;(* Based on your description, I made rc 30×10^-6 bigger. *)
s = NDSolve[{D[r*O2'[r], r] - con*r*((O2[r])/(O2[r] + km)) == 0, O2'[rc] == 0, O2[ro] == po},
            O2, {r, ro, rc}, Method -> {"Shooting", "StartingInitialConditions" -> 
                                                    {O2[rc] == 0, O2'[rc] == 0}}];      
{O2[rc], O2'[rc]} /. s
Plot[O2[t] /. s, {t, ro, rc}, PlotRange -> All]

{{0.167193, 0.}}

enter image description here

O2[rc] is a little different from what you give, I think it doesn't hurt and believe that it can be improved by more adjustment but I don't have time now. (According to my experience, equations which have to turn to shooting method are often sensitive. )

There're many posts about shooting method in this site, you can have a search.

  • $\begingroup$ That is excellent, thank you - shooting method works fine. As you do point out, it is extremely sensitive; at least this narrows the range I have to look in. Thanks again! $\endgroup$
    – DRG
    Oct 22, 2013 at 9:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.