# NDSolve::ndsz error: singularity or stiff system suspected | Coupled ODEs [closed]

I am new to Mathematica. Trying to solve the below coupled-ODEs (eqnA, eqnB, eqnC). However, getting error NDSolve::ndsz: At x == 0.8849324488967629, step size is effectively zero; singularity or stiff system suspected. Referring to other related posts, tried solutions like changing the method to StiffnessSwitching, and other options like WorkingPrecision, etc, but nothing is working.

g  = 5/3;
eqnA = (y2[x]-x)x y1'[x] + (x y2'[x]+ y2[x])y1[x] == 0;
eqnB= -3/2  y1[x]y2[x]  + (y2[x] - x)y1[x]y2'[x] + y3'[x] == 0;
eqnC = -3 y1[x] y3[x] + (y2[x] - x)(y1[x]y3'[x] - g y3[x] y1'[x] )== 0;
sol=NDSolve[{eqnA,eqnB, eqnC,y1[1]== (g+1)/(g-1),y2[1]== 2/(g+1), y3[1] == 2/(g+1)},{y1,y2,y3},{x,0,1}];
Plot[Evaluate[{y1[x],y2[x],y3[x]}/.sol[[1]]],{x,0,1}]


Expecting a plot somewhat like this.

• Solving for the derivatives shows that they become singular when 2 x (3 x^2 y1[x] - 6 x y1[x] y2[x] + 3 y1[x] y2[x]^2 - 5 y3[x])` vanishes. Check whether this is the cause of the integration stopping. Nov 19, 2021 at 21:06
• I just checked, and that is exactly what is happening. This is a math issue, not a Mathematica issue. Nov 19, 2021 at 21:15
• @bbgodfrey Ok! But how we can solve this problem? Nov 20, 2021 at 4:39
• @twilight Is there exact solution you shown? Nov 20, 2021 at 6:30
• The coefficients of the ODEs weren't right. When I solved them again, I got the right coefficients and the desired plot using them. Your comments have been more than helpful. Thank you very much. Dec 1, 2021 at 20:05