Singular ODEs with NDSolve [closed]

I have a system of singular ODEs:

eq1 =  y1'[x] y1[x] - x y1'[x]^2 + x y1''[x] y1[x] == c1 y2[x];
eq2 = y2'[x] y2[x] - x y2'[x]^2 + x y2''[x] y2[x] == c2 y1[x];

This system has a regular (at $x = 0$) analytical solution:

$y_1(x) = 1 + c_1 x + \frac{1}{4} c_1 c_2 x^2, \quad y_2(x) = 1 + c_2 x + \frac{1}{4} c_1 c_2 x^2$.

However, I am unable to use NDSolve to find this solution. I am trying to do the following:

c1 = 0.25;
c2 = 0.5;

NDSolve[{eq1 == 0, eq2 == 0, y1 ==  1, y2 == 1, y1' == c1, y2' == c2},
{y1, y2}, {x, 0, 1000}, Method -> {"EquationSimplification" -> "Residual"},
AccuracyGoal -> 10]

But I arrive at

NDSolve::icfail: Unable to find initial conditions that satisfy the residual function
within specified tolerances. Try giving initial conditions for both values and
derivatives of the functions. >>

I am trying to use "Residual" method because equations are singular at $x = 0$. I also do not understand why NDSolve cannot 'find' (?) initial conditions that satisfy the residual function. What does it mean? I am giving these conditions explicitly, but it asks me to give initial conditions for both values and derivatives which I exactly do.

I tried to avoid point $x = 0$ by setting some epsilon, say

eps = 0.00000001

and start integrating from this point (with initial conditions on eps provided), but it does not help.

Can anybody tell me how I could resolve this problem and integrate this system successfully using the NDSolve function?

closed as off-topic by Michael E2, LCarvalho, Coolwater, Sektor, m_goldbergDec 3 '17 at 22:32

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Michael E2, LCarvalho, Coolwater, Sektor, m_goldberg
If this question can be reworded to fit the rules in the help center, please edit the question.

• You have == in your equations twice. Once in your definitions, and then again when you set them equal to zero. – aardvark2012 Dec 3 '17 at 10:17
• With @aardvark2012's fix, everything seems to work fine. No errors, solution returned. – Michael E2 Dec 3 '17 at 17:13

You should use ParametricNDSolve! Because of the parameters c1, c2 N(umerical!)DSolve can't solve the problem:

sol = ParametricNDSolve[{eq1 , eq2 ,
y1 == 1, y2 == 1, y1' == c1, y2' == c2},
{y1, y2}, {x, 0, 1000}, {c1, c2},
Method -> {"EquationSimplification" -> "Residual"},AccuracyGoal -> 10]

With give c1 = 0.25;c2 = 0.5; NDSolve finds the solution directly:

erg = NDSolve[{eq1 , eq2 , y1 == 1, y2 == 1, y1' == c1,
y2' == c2}, {y1, y2}, {x, 0, 1000},
Method -> {"EquationSimplification" -> "Residual"},
AccuracyGoal -> 10]
Plot[{y1[x], y2[x]} /. erg[]  , {x, 0, 10}] without error. The initial conditions are as requested!

• Thank you for the answer. But I fixed the values c1 and c2. What might be wrong?.. – newt Dec 3 '17 at 17:03
• Nevertheless, if I use ParametricNDSolve and substitute the values of parameters that I used in the question, I get the same error: Unable to find initial conditions that satisfy the residual function \ within specified tolerances. Try giving initial conditions for both \ values and derivatives of the functions. – newt Dec 3 '17 at 17:09
• Clear[c1,c2] before using ParametricNDSolve. Then define their values before using the NDSolve part. Worked for me. – Bill Watts Dec 3 '17 at 21:41