# ND Solve not working : Step size zero, singularity/stiff system suspected

I am trying to solve a set of four coupled equations through NDSolve. However, I am getting the error "step size is effectively zero; singularity or stiff system suspected.." I think this is a common error when there exists very large or very big numbers, but I can't see any such numbers in my system. So, I am worried what could be the reason ? The set of equations are given below :

eqn1 = M'[a] +
Log[10]/Sqrt[1/3  (B[a] 10^(-3 a) + R [a] 10^(-4 a) )] (1/M[a]) ^2;
eqn2 = R'[a] + 10^a/M[a] M'[a] B[a];
eqn3 = B'[a] - 1/M[a] M'[a] B[a];
eqn4 =  s'[a] + 3 s[a] Log[10] +
1/  (s[a]^((1/3))) B[a]/(M[a] 10^a) M'[a];
sol1 = NDSolve[{eqn1 == 0 , eqn2 == 0, eqn3 == 0,
eqn4 == 0, M[2] == 2, R[2] == 1,
B[2] == 0.001, s[2] == 0.1}, {M, R, B , s}, {a, 2, 8}]

• Have you plotted the solution? If the solution has a singularity, it may be a feature of the system and there is nothing to do about it. (E.g. how would you get rid of the pole at $x=0$ of $y=1/x$?). The singularity may also be due to the system being unstable in the neighborhood of your IVP. Related: (39028), (232756). – Michael E2 Jan 12 at 16:45

Seems for all positive initial conditions (m2) for M, system has a singularity. But for negative m2, you get a full a-range solution.

ndsol[b2_, m2_, r2_, s2_] :=
First@NDSolve[{eqn1 == 0, eqn2 == 0, eqn3 == 0, eqn4 == 0,
B[2] == b2, M[2] == m2, R[2] == r2, s[2] == s2}, {B, M, R, s}, {a,
2, 8}]

Plot[Evaluate[M[a] /. ndsol[.001, 2, 1, .1]], {a, 2,
2 + .9999*(puf[[1, 2, 1, 1, 2]] - 2)}]

Manipulate[
Plot[Evaluate[M[a] /. (puf = ndsol[b2, m2, r2, s2] // Quiet)], {a, 2,
2 + .9999*(puf[[1, 2, 1, 1, 2]] - 2)}], {{b2, .001}, -3, 3,
Appearance -> "Labeled"}, {{m2, -2}, -5, 5,
Appearance -> "Labeled"}, {{r2, 1}, -5, 5,
Appearance -> "Labeled"}, {{s2, .001}, -1, 1,
Appearance -> "Labeled"}]