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.
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. $\endgroup$