# Solving a nonlinear system of 2nd order ODEs with NDSolve

I'm trying to solve a 2nd order related ODE system with five dependent variables.

The code is:

    Nc = 8; Nn = Nc-4; Pt = 0.495; R = 0.08206; Ru = 8.314; Cebo = 0.992; Cwo = 5.76; Pebo = 0.073;
Pwo = 0.422;  yebo = 0.1428; ywo = 0.8571; T0 = 893.15;
yeb[z_] := (Ceb[z] R T[z])/(Pt 1000); yw = ywo; yst[z_] :=  Abs[(1 - Ceb[z] - ywo)/2];
c2[z_] := (yw Pt 1000)/(R T[z]); c3[z_] := (yst[z] Pt 1000)/(R T[z]); c4[z_] := c3[z];
Mmeb = 106.17; Mmw = 18; Mmst = 104.15; Mmh2 = 2; Mmbe = 78.11; Mmet = 28.05; Mmto = 92.14; Mmme= 16.04;
M[1] = Mmeb; M[2] = Mmw; M[3] = Mmst; M[4] = Mmh2; M[5] = Mmbe; M[6] = Mmet; M[7] = Mmto; M[8] = Mmme;
j = 0;
Do[i = 1; j = j + 1;
M[i_,j_] := 2*(1/M[i] + 1/M[j])^-1, Nc];
Teb = 409.15; Tw = 373.15; Tst= 418.15; Th2 = 20.27; Tbe = 353.25; Tet = 169.5; Tto = 384.00; Tme= 111.65;
Tb[1] = Teb;  Tb[2] = Tw; Tb[3] = Tst; Tb[4] = Th2; Tb[5] = Tbe; Tb[6] = Tet; Tb[7] = Tto;Tb[8] = Tme;
Veb = 139.24; Vw = 18.789; Vst= 131.27; Vh2 = 28.16; Vbe = 96.017; Vet = 49.29; Vto = 118.29; Vme= 35.64;
Vm[1] = Veb; Vm[2] = Vw; Vm[3] = Vst; Vm[4] =Vh2; Vm[5] = Vbe; Vm[6] = Vet; Vm[7] = Vto; Vm[8] = Vme;
Deb = 0.58; Dw = 1.8546; Dst = 0.13; Dh2 = 0; Dbe = 0; Dett = 0; Dto = 0.36; Dme = 0;
Dp[1] = Deb; Dp[2] = Dw; Dp[3] = Dst; Dp[4] = Dh2; Dp[5] = Dbe; Dp[6] = Det; Dp[7] = Dto; Dp[8] = Dme;
i = 0;
Do[i = i + 1;
DD[i_] = (1940*Dp[i]^2)/(Vm[i]*Tb[i]); SS[i_] = ((1.585*Vm[i])/(1 + 1.3*DD[i]^2))^(1/3); EK[i_] = 1.18 (1 + 1.3*DD[i]^2)*Tb[i];
CO[i_,z_] := (1.16145*(T[z]/EK[i])^-0.14874) + 0.52487*(Exp[-0.77320 (T[z]/EK[i])]) + 2.16178*(Exp[-2.43787 (T[z]/EK[i])]);
Mu[i_,z_] := (2.6709 10^-6)*Sqrt[M[i]*T[z]]/((SS[i]^2)*CO[i,z]),Nc];
y[1,z_] := yeb[z]; y[2,z_] := yw; y[3,z_] := yst[z]; y[4,z_] := yst[z];
i = 0; j = 0;
Do[i = i + 1;
PHI[i,i,z] = 1;
Do[j = j + 1;
PHI[i_,j_,z_] := (1/Sqrt[8]*(1 + M[i]/M[j])^-0.5)*(1 + ((Mu[i,z]/Mu[j,z])^0.5)*(M[j]/M[i])^0.25)^2, Nc],Nc];
Mum[z_] := \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$4$$]$$(\*FractionBox[\(y[i, z]*Mu[i, z]$$, $$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$4$$]$$\(y[j,z]$$$$\$$$$PHI[i, j, z]$$$$\$$\)\)])\)\);
Daeb := 8*15.9 + 10*2.31; Daw := 4.54 + 2*2.31; Dast := 8*15.9 + 8*2.31; Dah2 := 2*2.31;
Da[1] = Daeb; Da[2] = Daw; Da[3] = Dast; Da[4] = Dah2;
i = 1; j = 0;
Do[j = j + 1;
Dc[i_,j_,z_] := (0.00143*T[z]^1.75)/(Pt M[i,j]*(Da[i]^(1/3) + Da[j]^(1/3))^2), Nc];
aa = 1; bb = 0; cc = 1; dd = 1;
ZZ= (cc + dd - aa - bb)/aa;
Dm[z_] := 1/(1/(1 + ZZ*y[1,z])*((1/Dc[1,2,z])*(y[2,z] - bb/aa y[1,z]) + ((1/Dc[1,3,z])*y[3,z] + cc/aa y[1,z]) + ((1/Dc[1,4,z])*(y[4,z] + dd/aa y[1,z]))));
n = 3;
Cp[1,z_] := (5 + n)*(1/2)*((8.314*T[z]^0)/M[1]); Cp[2,z_] := 143.05 - 183.54*(T[z]/100)^0.25 + 82.751*(T[z]/100)^0.5 - 3.6989*(T[z]/100)^1; Cp[3,z_] := (5 + n)*(1/2)*(8.314 /M[3])*T[z]^0 ; Cp[4,z_] := 56.505 - 702.24 (T[z]/100)^-0.75 + 1165.0 (T[z]/100)^-1 - 560.70 (T[z]/100)^-1.5;
Cpm[z_] := \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$4$$]\ $$Cp[i, z]\ y[i,z]$$\);
i = 0;
j = 0;
Do[i = i + 1;
k[i_,z_] := ((8.3127*10^-2)*Sqrt[(T[z]/M[i])])/((SS[i]^2)*CO[i,z]) + 1.32*(Cp[i,z] - (5/2)*8.314/M[i])*((2.6709*10^-6)*Sqrt[M[i]*T[z]])/((SS[i]^2)*CO[i,z]),Nn];
km[z_] := \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$4$$]$$(\*FractionBox[\(y[i, z]\ k[i, z]$$, $$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$4$$]$$\(y[j,z]$$$$\$$$$PHI[i, j, z]$$$$\$$\)\)])\)\);
CDeb = -102.22; CDst = -104.56; CDh2 = -117.95;
Ea[1] = 175.38; Ea[2] = 296.29; Ea[3] = 474.76; Ea[4] = 213.78;
Aeb = (1.014*10^-5)/0.9869; Ast = (2.678*10^-5)/0.9869; Ah2 = (4.519*10^-7)/0.9869;
A[1] = 5.594*10^9; A[2] = 1.060*10^15; A[3] = 1.246*10^26; A[4] = 8.024*10^10;
i = 0;
Do[i = i + 1; kc[i_,z_] := A[i]*Exp[-Ea[i]/(Ru/1000 T[z])],Nn];
keb[z_] := Aeb*Exp[-CDeb/(Ru/1000 T[z])]; kst[z_] := Ast*Exp[-CDst/(Ru/1000 T[z])]; kh2[z_] := Ah2*Exp[-CDh2/(Ru/1000 T[z])];
KEQ[z_] := (y[3,z]*Pt*y[4,z]*Pt)/(y[1,z]*Pt);
r1[z_] := (kc[1,z]*keb[z]*(y[1,z]*Pt - (y[3,z]*Pt*y[4,z]*Pt)/(KEQ[z]*Pt)))/(1 + keb[z]*y[1,z]*Pt + kh2[z]*y[4,z]*Pt + kst[z]*y[3,z]*Pt)^2;
r2[z_]:= (kc[2,z]*keb[z]*y[1,z]*Pt)/(1 +keb[z]*y[1,z]*Pt + kh2[z]*y[4,z]*Pt + kst[z]*y[3,z]*Pt)^2;
r3[z_] := (kc[3,z]*keb[z]*y[1,z]*Pt*kh2[z]*y[4,z]*Pt)/(1 + keb[z]*y[1,z]*Pt + kh2[z]*y[4,z]*Pt + kst[z]*y[3,z]*Pt)^2;
CDH1 = 124.83; CDH2=101.50; CDH3=-65.06;
u=((0.005 + 0.00585)/2)*3600; q = (35.52 + 31.8)/2; L=6;
eqn1 = Dm[z]*Ceb''[z] - u*Ceb'[z] == 1442*r1[z];
eqn2 = km[z]*T''[z] - q*Cpm[z]*T'[z] == 1442*(r1[z]*(-CDH1) + r2[z]*(-CDH2) + r3[z]*(-CDH3));
initConds={T[0] == 893.15, T'[L] == 0, Ceb[0] == 0.992, Ceb'[L] == 0 };
eqns=Join[{eqn1,eqn2}];
NDSolve[{eqns,initConds},{Ceb[z],T[z]},{z,0,6,0.01}]


I'm using to use NDSolve as shown below, but the code did not solve the system

Assuming that Mathematica can solve my system, how could I write code using NDSolve or any other method?

• C and D are predefined functions in Mathematica and can't be used as your own functions without grief. Perhaps rename those CapitalC and CapitalD. Next, NDSolve expects every variable and function you aren't solving for to have a predefined, preferably numeric value, so u,a,r1,b,r2,r3,d don't have values that the reader can use to try in his own worksheet and see if his potential answer actually works. Perhaps this will give you enough correct the first wave of things in your code. – Bill Dec 30 '18 at 19:13
• Thanks Bill! I would write the corresponding functions of r1,r2,r2 but the problem is that those functions are related to other in "waterfall" method. If so, I should write all codec here and I think it's not the way to ask. Those functions can be whatever only if they depends of T[z]. For the rest, "u" is equal to 0.005, "q" is equal to 32, "a"=1442, "b"=124, "c"=100 and "d"=-65. – atavistic mephit Dec 30 '18 at 19:32

I changed the code a bit, but the final system of equations does not differ from that obtained by the author's code. To solve the system we use the method of the false transient.

Nc = 8; Nn =
Nc - 4; Pt = 0.495; R = 0.08206; Ru = 8.314; Cebo = 0.992; Cwo = \
5.76; Pebo = 0.073;
Pwo = 0.422;  yebo = 0.1428; ywo = 0.8571; T0 = 893.15;
yeb[z_] := (Ceb[z] R T[z])/(Pt 1000); yw = ywo;
yst[z_] :=  Abs[(1 - Ceb[z] - ywo)/2];
c2[z_] := (yw Pt 1000)/(R T[z]); c3[z_] := (yst[z] Pt 1000)/(R T[z]);
c4[z_] := c3[z];
Mmeb = 106.17; Mmw = 18; Mmst = 104.15; Mmh2 = 2; Mmbe = 78.11; Mmet \
= 28.05; Mmto = 92.14; Mmme = 16.04;
M[1] = Mmeb; M[2] = Mmw; M[3] = Mmst; M[4] = Mmh2; M[5] = Mmbe;
M[6] = Mmet; M[7] = Mmto; M[8] = Mmme;

Table[
M1[i, j] = 2*(1/M[i] + 1/M[j])^-1, {i, 1, Nc}, {j, 1, Nc}];
Teb = 409.15; Tw = 373.15; Tst = 418.15; Th2 = 20.27; Tbe = 353.25; \
Tet = 169.5; Tto = 384.00; Tme = 111.65;
Tb[1] = Teb;  Tb[2] = Tw; Tb[3] = Tst; Tb[4] = Th2; Tb[5] = Tbe;
Tb[6] = Tet; Tb[7] = Tto; Tb[8] = Tme;
Veb = 139.24; Vw = 18.789; Vst = 131.27; Vh2 = 28.16; Vbe = 96.017; \
Vet = 49.29; Vto = 118.29; Vme = 35.64;
Vm[1] = Veb; Vm[2] = Vw; Vm[3] = Vst; Vm[4] = Vh2; Vm[5] = Vbe;
Vm[6] = Vet; Vm[7] = Vto; Vm[8] = Vme;
Deb = 0.58; Dw = 1.8546; Dst = 0.13; Dh2 = 0; Dbe = 0; Dett = 0; Dto \
= 0.36; Dme = 0;
Dp[1] = Deb; Dp[2] = Dw; Dp[3] = Dst; Dp[4] = Dh2; Dp[5] = Dbe;
Dp[6] = Det; Dp[7] = Dto; Dp[8] = Dme;

Do[
DD[i] = (1940*Dp[i]^2)/(Vm[i]*Tb[i]);
SS[i] = ((1.585*Vm[i])/(1 + 1.3*DD[i]^2))^(1/3);
EK[i] = 1.18 (1 + 1.3*DD[i]^2)*Tb[i];
CO[i][z_] := (1.16145*(T[z]/EK[i])^-0.14874) +
0.52487*(Exp[-0.77320 (T[z]/EK[i])]) +
2.16178*(Exp[-2.43787 (T[z]/EK[i])]);
Mu[i][z_] := (2.6709 10^-6)*
Sqrt[M[i]*T[z]]/((SS[i]^2)*CO[i][z]);, {i, 1, Nc}];
y[1, z_] := yeb[z]; y[2, z_] := yw; y[3, z_] := yst[z];
y[4, z_] := yst[z];

Do[
PHI[i, i][z_] := 1;
Do[
PHI[i, j][
z_] := (1/
Sqrt[8]*(1 +
M[i]/M[j])^-0.5)*(1 + ((Mu[i][z]/Mu[j][z])^0.5)*(M[j]/
M[i])^0.25)^2;, {i, 1, Nc}], {j, 1, Nc}];
Mum[z_] := \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$4$$]$$( \*FractionBox[\(y[i, z]*\(Mu[i]$$[z]\), $$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$4$$]$$\(y[j, z]$$$$\$$$$\(PHI[i, j]$$[z]\)$$\$$\)\)])\)\);
Daeb := 8*15.9 + 10*2.31; Daw := 4.54 + 2*2.31; Dast :=
8*15.9 + 8*2.31; Dah2 := 2*2.31;
Da[1] = Daeb; Da[2] = Daw; Da[3] = Dast; Da[4] = Dah2;

Table[
Dc[j] = (0.00143*T[z]^1.75)/(Pt *
M1[1, j]*(Da[1]^(1/3) + Da[j]^(1/3))^2), {j, 1, Nc}];
aa = 1; bb = 0; cc = 1; dd = 1;
ZZ = (cc + dd - aa - bb)/aa;
Dm[z_] := 1/(1/(1 +
ZZ*y[1, z])*((1/Dc[2])*(y[2, z] - bb/aa y[1, z]) + ((1/Dc[3])*
y[3, z] +
cc/aa y[1, z]) + ((1/Dc[4])*(y[4, z] + dd/aa y[1, z]))));
n = 3;
Cp[1, z_] := (5 + n)*(1/2)*((8.314*T[z]^0)/M[1]);
Cp[2, z_] :=
143.05 - 183.54*(T[z]/100)^0.25 + 82.751*(T[z]/100)^0.5 -
3.6989*(T[z]/100)^1;
Cp[3, z_] := (5 + n)*(1/2)*(8.314 /M[3])*T[z]^0 ;
Cp[4, z_] :=
56.505 - 702.24 (T[z]/100)^-0.75 + 1165.0 (T[z]/100)^-1 -
560.70 (T[z]/100)^-1.5;
Cpm[z_] := \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$4$$]\ $$Cp[i, z]\ y[i, z]$$\);
i = 0;
j = 0;
Do[
k[i][z_] := ((8.3127*10^-2)*Sqrt[(T[z]/M[i])])/((SS[i]^2)*
CO[i][z]) +
1.32*(Cp[i, z] - (5/2)*8.314/M[i])*((2.6709*10^-6)*
Sqrt[M[i]*T[z]])/((SS[i]^2)*CO[i][z]), {i, 1, Nn}];
km[z_] := \!$$\*UnderoverscriptBox[\(\[Sum]$$, $$i = 1$$, $$4$$]$$( \*FractionBox[\(y[i, z]\ \(k[i]$$[z]\), $$\*UnderoverscriptBox[\(\[Sum]$$, $$j = 1$$, $$4$$]$$\(y[j, z]$$$$\$$$$\(PHI[i, j]$$[z]\)$$\$$\)\)])\)\);
CDeb = -102.22; CDst = -104.56; CDh2 = -117.95;
Ea[1] = 175.38; Ea[2] = 296.29; Ea[3] = 474.76; Ea[4] = 213.78;
Aeb = (1.014*10^-5)/0.9869; Ast = (2.678*10^-5)/
0.9869; Ah2 = (4.519*10^-7)/0.9869;
A[1] = 5.594*10^9; A[2] = 1.060*10^15; A[3] = 1.246*10^26;
A[4] = 8.024*10^10;

Table[ kc[i] = A[i]*Exp[-Ea[i]/(Ru/1000 T[z])], {i, 1, Nn}];
keb[z_] := Aeb*Exp[-CDeb/(Ru/1000 T[z])];
kst[z_] := Ast*Exp[-CDst/(Ru/1000 T[z])];
kh2[z_] := Ah2*Exp[-CDh2/(Ru/1000 T[z])];
KEQ[z_] := (y[3, z]*Pt*y[4, z]*Pt)/(y[1, z]*Pt);
r1[z_] := (kc[1]*
keb[z]*(y[1, z]*Pt - (y[3, z]*Pt*y[4, z]*Pt)/(KEQ[z]*Pt)))/(1 +
keb[z]*y[1, z]*Pt + kh2[z]*y[4, z]*Pt + kst[z]*y[3, z]*Pt)^2;
r2[z_] := (kc[2]*keb[z]*y[1, z]*
Pt)/(1 + keb[z]*y[1, z]*Pt + kh2[z]*y[4, z]*Pt +
kst[z]*y[3, z]*Pt)^2;
r3[z_] := (kc[3]*keb[z]*y[1, z]*Pt*kh2[z]*y[4, z]*
Pt)/(1 + keb[z]*y[1, z]*Pt + kh2[z]*y[4, z]*Pt +
kst[z]*y[3, z]*Pt)^2;
CDH1 = 124.83; CDH2 = 101.50; CDH3 = -65.06;
u = ((0.005 + 0.00585)/2)*3600; q = (35.52 + 31.8)/2; L = 6; t0 = 2;
eqns = { Dm[z]*Ceb''[z] - u*Ceb'[z] - 1442*r1[z],
km[z]*T''[z] - q*Cpm[z]*T'[z] -
1442*(r1[z]*(-CDH1) + r2[z]*(-CDH2) + r3[z]*(-CDH3))} /. {T[z] ->
U[t, z], Ceb[z] -> W[t, z], T'[z] -> D[U[t, z], z],
T''[z] -> D[U[t, z], z, z], Ceb'[z] -> D[W[t, z], z],
Ceb''[z] -> D[W[t, z], z, z]};
bc = {U[t, 0] == 893.15, Derivative[0, 1][U][t, L] == 0,
W[t, 0] == 0.992, Derivative[0, 1][W][t, L] == 0 };
ic = {U[0, z] == 893.15, W[0, z] == 0.992}; L = 6; t0 = 2;
sol = NDSolve[{eqns == {D[W[t, z], t], D[U[t, z], t]}, bc, ic}, {U,
W}, {t, 0, t0}, {z, 0, L},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 80, "MaxPoints" -> 100,
"DifferenceOrder" -> 2}}, MaxSteps -> 10^6];
{Plot3D[Evaluate[U[t, z] /. sol], {t, 0, t0}, {z, 0, L}],
Plot3D[Evaluate[W[t, z] /. sol], {t, 0, t0}, {z, 0, L}],
Plot[Evaluate[U[t0, z] /. sol], {z, 0, L}, AxesLabel -> {"z", "T"}],
Plot[Evaluate[W[t0, z] /. sol], {z, 0, L},
AxesLabel -> {"z", "Ceb"}],
Plot[Evaluate[U[t0, z] /. sol], {z, 0.95*L, L},
AxesLabel -> {"z", "T"}],
Plot[Evaluate[W[t0, z] /. sol], {z, 0.95*L, L},
AxesLabel -> {"z", "Ceb"}]}


Compare this character by character with what you have. There is a missing ] in your code that I guessed where to insert, there appear to be differences in capitalization and naming, etc. Maybe this will help you get started. Please run this code as is in a fresh empty notebook and verify if it works for you. Then make small changes, one at a time, trying to bring this closer to your actual problem and verify that it still works after each tiny change, finding and fixing the source of any errors that it displays before making the next tiny change. If you can provide corrections then I will try to edit those into this

L=9;u=0.005;a=1442;b=124;c=100;d=65;p=5;q=32;
CapitalD[z_]:=9;r2[z_]:=11;r3[z_]:=13;
k[z_]:=17;kc1[z_]:=19;keb[z_]:=23;
keq[z_]:=29;kh2[z_]:=31;kst[z_]:=37;
y1[z_]:=41;y3[z_]:=43;y4[z_]:=47;
r1[z_]:=(kc1[z]*keb[z]*(y1[z]*p-(y3[z]*p+y4[z]*p)/keq[z]*p))/
(1+keb[z]*y1[z*p]+kh2[z]*y4[z]*p+kst[z]*y3[z]*p)^2;
eqn1 = CapitalD[z]*CapitalC''[z]-u*CapitalC'[z]==a*r1[z];
eqn2 = k[z]*T''[z]-q*T'[z]==a*(r1[z]*b+r2[z]*c+r3[z]*CapitalD[z]);
initConds={T[0]==893.15, T'[L]==0,CapitalC[0]==0.992,CapitalC'[L]==0};
s = NDSolve[{eqn1, eqn2, initConds}, {T[z], CapitalC[z]}, {z, 0, 9}][[1]];
Plot[Evaluate[CapitalC[z]/.s],{z,0,9}]
Plot[Evaluate[T[z]/.s],{z,0,9}]

• Okay, will rewrite and post the code as soon as possible. – atavistic mephit Dec 31 '18 at 1:33
• Bill! I used your code in my worksheet and well, it plot something. At the moment this plot has no sense, but at least it plot something. By the way, in your code are lot of constants that should be variables as r2[z], r3[z], y1[z] etc.. – atavistic mephit Dec 31 '18 at 9:45
• Okay, using your code changing all constants to his respective functions, the program give me the next output: NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}. – atavistic mephit Dec 31 '18 at 10:42
• as NDSolve is not able to solve this system, Do u know any other method capable of it? – atavistic mephit Dec 31 '18 at 15:16
• okay, srry then! Didn't want to bother you. Thanks for the help anyways! – atavistic mephit Dec 31 '18 at 17:41

Assuming that previously all the needed functions and constant are defined, you can try

kc1[z_] := 1 + z
keb[z_] := 1 + z
y1[z_] := 1 + z
y3[z_] := 1 + z
y4[z_] := 1 + z
keq[z_] := 1 + z
kh2[z_] := 1 + z
kst[z_] := 1 + z
d[z_] := 1 + z
k[z_] := 1 + z
r1[z_] := (kc1[z]*
keb[z]*(y1[z]*p - (y3[z]*p + y4[z]*p)/keq[z]*p))/(1 +
keb[z]*y1[z]*p + kh2[z]*y4[z]*p + kst[z]*y3[z]*p)^2;
r2[z_] := r1[z]
r3[z_] := r1[z]
eqn1 = d[z]*CC''[z] - u*CC'[z] == a*r1[z];
eqn2 = k[z]*T''[z] - q*T'[z] == a*(r1[z]*b + r2[z]*c + r3[z]*d0);
initConds = {T[0] == 893.15, T'[L] == 0, CC[0] == 0.992, CC'[L] == 0};
parms = {p -> 1, q -> 1, a -> 1, b -> 1, c -> 1, d0 -> 1, u -> 1};
L = 10;
zmax = 9;
delz = 0.01;
eqns = {eqn1, eqn2} /. parms;
DE = Join[eqns, initConds];
s = NDSolve[DE, {T, CC}, {z, 0, zmax, delz}]
Plot[Evaluate[{T[z], CC[z]} /. s], {z, 0, zmax}]

• Hello Cesareo, I also used your code before and the result is the same. NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}. – atavistic mephit Dec 31 '18 at 10:49
• @atavisticmephit To avoid misunderstandings in this edit, I defined the declared functions and associated values to the parameters used. I hope this helps. – Cesareo Dec 31 '18 at 11:21
• @atavisticmephit This code as it is, runs without problems at my computer with MATHEMATICA v 11.0 – Cesareo Dec 31 '18 at 11:27
• I going to try it, let's see if it works but I don't understant those parameters. Why all of them are 1? – atavistic mephit Dec 31 '18 at 11:28
• Well, using it i have a solution. The program draws a plot, with no sense, but a plot. There are things that I don't understant as {p -> 1, 33.66 -> 1, a -> 1, b -> 1, c -> 1, d0 -> 1, 19.53 -> 1} ( output message). – atavistic mephit Dec 31 '18 at 11:35