I am trying to solve the following differential equation

EOM=(-Piecewise[{{1.157793814828612*^106, r >= 0 && r <= 1.0408060341965232*^-18}, {(9.734399999999999*^36*(2.7912782245116803*^-35 - 3*r^2))/
         (2.7912782245116803*^-35 + r^2)^3, r > 1.0408060341965232*^-18}}, 0])*f[r] + (7*Derivative[1][f][r])/r + Derivative[2][f][r] == 0

with initial conditions

f[0] == 1
Derivative[1][f][0] == 0

numerically using NDSolve

NDSolve[{EOM, {f[0] == 1, Derivative[1][f][0] == 0}}, {f[r]}, {r, 0, 100}]

But I get

Power::infy: Infinite expression 1/0 encountered.

Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered.

NDSolve::ndnum: Encountered non-numerical value for a derivative at r == 0.

I searched other related posts but could not solve the problem. I will be thankful if you help me with this.

  • $\begingroup$ Your EOM contains a term (7*f'[r])/r which results in Infinite expression 1/0 encountered at the initial conditions (r = 0). $\endgroup$
    – Bob Hanlon
    Commented Jun 24, 2021 at 14:00
  • $\begingroup$ Thank you for your reply, I actually tried multiplying the whole equation by r but that did not help and I still get the same error. @Bob Hanlon $\endgroup$
    – amy
    Commented Jun 24, 2021 at 15:14
  • $\begingroup$ A little besides the point, but the numbers you have seem a bit wacky - they span 141 orders of magnitude which will likely create some numerical issues for you. Are you sure they're correct? $\endgroup$
    – bRost03
    Commented Jun 24, 2021 at 16:41
  • $\begingroup$ Change the initial conditions to {f[epsilon] == 1, f'[epsilon] == 0} for a sufficiently small epsilon $\endgroup$
    – Bob Hanlon
    Commented Jun 24, 2021 at 17:18
  • $\begingroup$ Yes, they are correct. The function pieces are of the same order of magnitude and also all the terms added up are of the same order of magnitude. So, in principle the equation is Okay. @bRost03 $\endgroup$
    – amy
    Commented Jun 24, 2021 at 17:23

2 Answers 2


This may be an answer to an approximation to your question. I am making assumptions.

In addition to the problem at r=0, the derivative of your piecewise function is zero when r=0. The this gives you a solution:

DSolve[{(7 Derivative[1][f][r])/r + (f^\[Prime]\[Prime])[r] == 0, 
  f'[0] == 0, f[0] == 1}, f, r]

and so the integrator isn't going to behave well near zero.

It look to me like you are approximating a potential that looks like this:

pot = a/r^5 - b/r^3

which gives you a force that looks like this.

force = -Simplify[D[pot, r]]

which looks like your piecewise function. Please forgive me for making assumptions about what you are doing.


This problem is begging to be rescaled to non-dimensional units. Let's suppose the mininum potential (-eMin) is located at rMin

solNorm = 
 Solve[{rMin == SolveValues[force == 0, r][[2]], 
    eMin == pot /. r -> rMin}, {a, b}][[1]]
    eNorm = Simplify[((pot /. r -> rho rMin)/eMin) /. solNorm]

That gives a normalized force term (to get units, multiply by Emin/rMin)

fNorm = - D[eNorm, rho]

With the assumption, your equation of motion becomes

eom = phi[rho] fNorm + 7 D[phi[rho], rho]/rho + D[phi[rho], {rho, 2}]

Note that starting the equation of motion at pho=0 is going to make phi fly off to infinity. Using @Bob Hanlon's advice, here is a function to get a numerical solution for different starting points.

phiSol[sp_] := 
 phiSol[sp] = 
  NDSolveValue[{eom == 0, phi'[sp] == 0, phi[sp] == 1}, 
   phi, {rho, sp, 2}]

For example:


We will need to extract the domain of the InterpolationFunction


The behavior can be investigated for different initial points:

 With[{interpolator = phiSol[sp], 
   domain = Flatten[InterpolatingFunctionDomain[phiSol[sp]]]},
  Plot[interpolator[t], {t, domain[[1]], domain[[2]]}, Frame -> True, 
   FrameLabel -> {"r/rMin", "phi"}]
 {{sp, .5}, .01, 1}

example Note that the length unit is rMin which is very small in your original problem.


I encountered similar problem. So I use answer space to tell you my numerical experiments.

Suppose that the solution $f(r)$ is even in $r$, then we may define $F(r^2):=f(r)$. Thus $$\partial_{rr}f(r)=\partial_{rr}F(r^2)=2F'(q)+4qF''(q),\quad q=r^2\tag{1}$$ $$\partial_{r}f(r)=\partial_{r}F(r^2)=2rF'(q),\quad\tfrac1r\partial_{r}f(r)=\tfrac1r\partial_{r}F(r^2)=2F'(q),\quad q=r^2\tag{2}$$

Thus EOM becomes a second order ODE for $F(q)$: $$0=A(q)F(q)+7*2F'(q)+2F'(q)+4qF''(q),\quad q=r^2\tag{3}$$

When I use the boundary condition $F(0)=1,F'(0)=0$ to NDSolve (3), I also got "Infinity expression 1/0 encountered". This is probably due to the term $4qF''(q)$. I then used boundary condition, suggested by
Bob Hanlon, $F(\epsilon)=1,F'(\epsilon)=0$ and was able to get solution but with different error message "At q == 1.0832772008198944`*^-38, step size is effectively zero; singularity or stiff system suspected". In theory the boundary condition $\partial_{r}f(0)==2rF'(q)|_{q=r^2=0}$ is automatically satisfied, no additional requirement is necessary on $F'(0)$. But if I do not specify $F'(0)$, the NDSolve won't work.

Another thing you may try is to substitute $F(q)=a_0 +q G(q)$ to get an ODE for $G(q)$.

Hope it helps -Mike


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.