# Numerically solving coupled 2rd ODE with tiny numbers

I am solving the initial value problem of two coupled ODEs numerically. The equations are

eqs = {(2.74315901125834 - 2.771065822236567/x^2) y1[x] + 2.743156268102072*^-6 y2[x] +
2.7989379082285724*^-18 y1''[x] == 0,
1. y1[x] + 1.0000000000634328 y2[x] + 1.0203348386583513*^-12 y2''[x] == 0,
y1[0.5] == 0, y2[0.5] == 1, y1'[0.5] == 0, y2'[0.5] == 0}


and I want to solve it with NDSolve

sol = NDSolve[eqs, {y1, y2}, {x, 0.5, 1.5}]


However, the results returned with errors

NDSolve::ndsz: At x == 0.5000000000000003, step size is effectively zero; singularity or stiff system suspected.


and the result is quite large, which is wrong.

I think the problem is from the tiny number, with the orders of $$10^{-12}$$ and $$10^{-18}$$.

I looked for solutions on the Internet and learned that using "ImplicitRungeKutta" method may make sense. So I tried

sol = NDSolve[eqs, {y1, y2}, {x, 0.5, 1.5},
Method -> {"FixedStep",
Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10}},
StartingStepSize -> 1/100]


and I tried different settings of StartingStepSize and DifferenceOrder, but they all failed. For large StartingStepSize such as $$10^{-3}$$, it returns

NDSolve::ndcf: Repeated convergence test failure at x == 1.004; unable to continue.


the same with large DifferenceOrder. If I choose "DifferenceOrder" -> 10" and "StartingStepSize -> 1/100", there are no errors, but the result is also quite large, which is in the order of $$10^{102}$$.

I wonder whether there are some methods to deal with tiny numbers in solving ODEs, or maybe I should use Python or Matlab and write the code using Runge-Kutta method which may perform better.

• you can try adding Method -> {"StiffnessSwitching"} and see what happens. But the correct way to handle this is to scale/normalize your ode itself. Commented May 8 at 4:25
• @Nasser Thank you for your comment. I have tried using Method -> {"StiffnessSwitching"}, but two hours have passed and the results haven't been calculated. I wonder if I can scale these ODEs because the tiny coefficient is not an overall coefficient.
• Eliminating the boundary conditions and replacing x^2 by c allows DSolve to give a very rough approximation involving terms like E^(9.98524*10^7 x). So, you should expect enormous answers from NDSolve`. Commented May 8 at 17:34