Dear all.

I am trying to use the ndsolve function to solve a system of ODEs, but I can't since I always obtain "singularity/stiffness system" errors. I would be very glad if someone can advice me or give me some assistance about this issue.

For more detail, I attach the Mathematica code and the obtained errors are summarized next:

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

InterpolatingFunction::dmval : Input value {0.8} lies outside the range of data in the interpolating function. Extrapolation will be used. InterpolatingFunction::dmval : Input value {0.9} lies outside the range of data in the interpolating function. Extrapolation will be used. InterpolatingFunction::dmval : Input value {1.} lies outside the range of data in the interpolating function. Extrapolation will be used.

General::stop : Further output of InterpolatingFunction::dmval will be suppressed during this calculation.

Thanks a lot. Best regards.

Manuel.

Mathematica code:

raftot = 3/1000;(*.003;*)
prmt5 = 3/1000;(*.003;*) (*0.0*)

e2tot = 3/10000;(*.0003;*)

mektot = 12/10;(*1.2;*)
mekpasetot = 3/10000;(*.0003;*)

erktot = 12/10;(*1.2;*)
erkpasetot = 12/100;(*.12;*)

a1 = 958.64;(*1000*)
a2 = 958.64;(*1000*)
a3 = 958.64;(*1000*)

a4 = 958.64;(*1000*)
a5 = 958.64;(*1000*)
a6 = 958.64;(*1000*)

a7 = 958.64;(*1000*)
a8 = 958.64;(*1000*)
a9 = 958.64;(*1000*)

a10 = 958.64;(*1000*)
aprmt5 = 482.94;

d1 = 150;
d2 = 150;
d3 = 150;
d4 = 150;
d5 = 150;
d6 = 150;
d7 = 150;
d8 = 150;
d9 = 150;
d10 = 150;

k1 = 150;
k2 = 150;
k3 = 150;
k4 = 150;
k5 = 150;
k6 = 150;
k7 = 150;
k8 = 150;
k9 = 150;
k10 = 150;

i = -1.5;
(*For[i=-6,i\[LessEqual]-6,i=i+diff,*)
e1tot = 10^(i);
alphaprmt5 = e1tot/10^(-1); (*alphaprmt5=((e1tot-mose1[t])/10^(-1)*)

s = NDSolve[{
    raf'[t] == -a1*raf[t]*e1[t] + d1*rafe1[t] + k2*rafpe2[t],
    rafe1'[t] == a1*raf[t]*e1[t] - (d1 + k1) rafe1[t],
    rafp'[t] == -a2*rafp[t]*e2[t] + d2*rafpe2[t] + 
      k1*rafe1[t] + (k3 + d3)*mekrafp[t] - 
      a3*rafp[t]*mek[t] + (k5 + d5)*mekprafp[t] - 
      a5*mekp[t]*rafp[t],(*-aprmt5*alphaprmt5*prmt5*mosstar[t],*)

    rafpe2'[t] == a2*rafp[t]*e2[t] - (d2 + k2)*rafpe2[t],
    mek'[t] == -a3*mek[t]*rafp[t] + d3*mekrafp[t] + 
      k4*mekpmekpase[t],
    mekrafp'[t] == a3*mek[t]*rafp[t] - (d3 + k3)*mekrafp[t],
    mekp'[t] == -a4*mekp[t]*mekpase[t] + d4*mekpmekpase[t] + 
      k3*mekrafp[t] + k6*mekppmekpase[t] + d5*mekprafp[t] - 
      a5*mekp[t]*rafp[t],
    mekpmekpase'[t] == 
     a4*mekp[t]*mekpase[t] - (d4 + k4)*mekpmekpase[t],
    mekprafp'[t] == a5*mekp[t]*rafp[t] - (d5 + k5)*mekprafp[t],
    mekpp'[t] == 
     k5*mekprafp[t] - a6*mekpp[t]*mekpase[t] + d6*mekppmekpase[t] - 
      a7*mekpp[t]*erk[t] + (d7 + k7)*erkmekpp[t] + (d9 + k9)*
       erkpmekpp[t] - a9*erkp[t]*mekpp[t],
    mekppmekpase'[t] == 
     a6*mekpp[t]*mekpase[t] - (d6 + k6)*mekppmekpase[t],
    erk'[t] == -a7*erk[t]*mekpp[t] + d7*erkmekpp[t] + 
      k8*erkperkpase[t],
    erkmekpp'[t] == a7*erk[t]*mekpp[t] - (d7 + k7)*erkmekpp[t],
    erkp'[t] == 
     k7*erkmekpp[t] - a8*erkp[t]*erkpase[t] + d8*erkperkpase[t] - 
      a9*erkp[t]*mekpp[t] + d9*erkpmekpp[t] + k10*erkpperkpase[t],
    erkpmekpp'[t] == a9*erkp[t]*mekpp[t] - (d9 + k9)*erkpmekpp[t],
    erkpp'[t] == -a10*erkpp[t]*erkpase[t] + d10*erkpperkpase[t] + 
      k9*erkpmekpp[t],
    erkperkpase'[t] == 
     a8*erkp[t]*erkpase[t] - (d8 + k8)*erkperkpase[t],
    erkpperkpase'[t] == 
     a10*erkpp[t]*erkpase[t] - (d10 + k10)*erkpperkpase[t],
    e1'[t] == -a1*raf[t]*e1[t] + (d1 + k1)*rafe1[t],
    e2'[t] == -a2*rafp[t]*e2[t] + (d2 + k2)*rafpe2[t],
    mekpase'[
      t] == -a4*mekp[t]*mekpase[t] + (d4 + k4)*mekpmekpase'[t] - 
      a6*mekpp[t]*mekpase[t] + (d6 + k6)*mekppmekpase[t],
    erkpase'[
      t] == -a8*erkp[t]*erkpase[t] + (d8 + k8)*erkperkpase'[t] - 
      a10*erkpp[t]*erkpase[t] + (d10 + k10)*erkpperkpase[t],

    raf[0] == raftot, rafe1[0] == 0, rafp[0] == 0, rafpe2[0] == 0, 
    mek[0] == mektot, mekrafp[0] == 0, mekp[0] == 0, 
    mekpmekpase[0] == 0, mekprafp[0] == 0, mekpp[0] == 0, 
    mekppmekpase[0] == 0, mekpase[0] == mekpasetot, erk[0] == erktot, 
    erkmekpp[0] == 0, erkp[0] == 0, erkpmekpp[0] == 0, erkpp[0] == 0, 
    erkperkpase[0] == 0, erkpperkpase[0] == 0, 
    erkpase[0] == erkpasetot, e1[0] == e1tot, e2[0] == e2tot},

   {raf, rafe1, rafp, rafpe2, mek, mekrafp, mekp, mekprafp, 
    mekpmekpase, mekpp, mekppmekpase, mekpase, erk, erkmekpp, erkp, 
    erkpmekpp, erkpp, erkperkpase, erkpperkpase, erkpase, e1, e2},

   {t, 0, 500}, MaxSteps -> 10000, 
   Method -> {StiffnessSwitching}]; 


Plot[raf[t] /. s, {t, 0, 0.8}, PlotRange -> All]

Plot[rafp[t] /. s, {t, 0, 0.8}, PlotRange -> All]

Plot[mekpp[t] /. s, {t, 0, 0.8}, PlotRange -> All]

Plot[erkpp[t] /. s, {t, 0, 0.8}, PlotRange -> All]