Hi all,
The DifferentialEquations’ suggested integration methods when trying to solve this stiff problem all yield the same error. The native Julia ones, that is. Sundial’s CVODE_BDF as well as lsoda are able to solve this with no problems.
My code:
using DifferentialEquations
using Plots
gr()
using LSODA
using Sundials
x0 = zeros(Float64, 35)
x0[1] = 5
x0[2] = 1
x0[3] = 5
x0[4] = 5
x0[5] = 1
x0[6] = 1
x0[7] = 1
x0[8] = 1
x0[9] = 1
x0[10] = 1
x0[11] = 1
x0[12] = 1
x0[13] = 1
x0[14] = 1
x0[15] = 1
x0[16] = 1
x0[17] = 1
x0[18] = 1
x0[19] = 1
x0[20] = 1
x0[21] = 1
x0[22] = 1
x0[23] = 1
x0[24] = 1
x0[25] = 1
x0[26] = 1
x0[27] = 1
x0[28] = 1
x0[29] = 1
x0[30] = 1
x0[31] = 1
x0[32] = 1
x0[33] = 1
x0[34] = 1
x0[35] = 6
p = zeros(Float64, 117)
p[1] = 1000
p[2] = 1000
p[3] = 1000
p[4] = 1000
p[5] = 1000
p[6] = 1000
p[7] = 1000
p[8] = 1000
p[9] = 1
p[10] = 1
p[11] = 0.7
p[12] = 0.7
p[13] = 0.7
p[14] = 0.7
p[15] = 0.2
p[16] = 0.2
p[17] = 1
p[18] = 1
p[19] = 0.7
p[20] = 0.7
p[21] = 0.05
p[22] = 1
p[23] = 0.5
p[24] = 0.7
p[25] = 0.5
p[26] = 2
p[27] = 0.7
p[28] = 0.7
p[29] = 2
p[30] = 0.2
p[31] = 0.2
p[32] = 1
p[33] = 0.7
p[34] = 0.2
p[35] = 0.21
p[36] = 1
p[37] = 0.7
p[38] = 1
p[39] = 0.2
p[40] = 0.2
p[41] = 0.05
p[42] = 0.01
p[43] = 0.7
p[44] = 1
p[45] = 0.5
p[46] = 1
p[47] = 0.04
p[48] = 0.01
p[49] = 0.01
p[50] = 0.01
p[51] = 1
p[52] = 1
p[53] = 0.5
p[54] = 1
p[55] = 1
p[56] = 0.5
p[57] = 0.1
p[58] = 0.1
p[59] = 0.5
p[60] = 1
p[61] = 1
p[62] = 1
p[63] = 0.3
p[64] = 1
p[65] = 1
p[66] = 0.7
p[67] = 0.4
p[68] = 0.7
p[69] = 1
p[70] = 0.7
p[71] = 0.7
p[72] = 2
p[73] = 0.2
p[74] = 2
p[75] = 1
p[76] = 0.7
p[77] = 0.7
p[78] = 0.2
p[79] = 0.7
p[80] = 0.5
p[81] = 0.21
p[82] = 0.2
p[83] = 0.2
p[84] = 0.01
p[85] = 2
p[86] = 2
p[87] = 2
p[88] = 2
p[89] = 2
p[90] = 2
p[91] = 2
p[92] = 2
p[93] = 1
p[94] = 1
p[95] = 1
p[96] = 1
p[97] = 1
p[98] = 0.5
p[99] = 0.7
p[100] = 0.7
p[101] = 2
p[102] = 0.2
p[103] = 2
p[104] = 1
p[105] = 1
p[106] = 0.7
p[107] = 1
p[108] = 0.7
p[109] = 1
p[110] = 0.7
p[111] = 1
p[112] = 0.7
p[113] = 2
p[114] = 2
p[115] = 0.7
p[116] = 0.2
p[117] = 0.2
function fun_b4(dx_dt, x, par, t)
p = zeros(347)
p[118] = 1
p[119] = 1
p[120] = 1
p[121] = 1
p[122] = 1
p[123] = 1
p[124] = 1
p[125] = 1
p[126] = 1
p[127] = 1
p[128] = 1
p[129] = 1
p[130] = 1
p[131] = 1
p[132] = 1
p[133] = 1
p[134] = 1
p[135] = 1
p[136] = 1
p[137] = 1
p[138] = 1
p[139] = 1
p[140] = 1
p[141] = 1
p[142] = 1
p[143] = 1
p[144] = 1
p[145] = 1
p[146] = 1
p[147] = 1
p[148] = 1
p[149] = 1
p[150] = 0
p[151] = 141.471
p[152] = 0.9
p[153] = 0.1
p[154] = 83402
p[155] = 127877
p[156] = 603.414
p[157] = 603.414
p[158] = 5.02845
p[159] = 21722.9
p[160] = 603.414
p[161] = 361200
p[162] = 673524
p[163] = 1.62679e+06
p[164] = 164390
p[165] = 36512.9
p[166] = 166804
p[167] = 39530
p[168] = 1810.24
p[169] = 127877
p[170] = 166201
p[171] = 1.94454e+06
p[172] = 5.02845
p[173] = 317755
p[174] = 0
p[175] = 532531
p[176] = 38926.6
p[177] = 40133.4
p[178] = 689982
p[179] = 4.45873e+06
p[180] = 2413.66
p[181] = 83402
p[182] = 603.414
p[183] = 603.414
p[184] = 2413.66
p[185] = 533135
p[186] = 100000
p[187] = 1
p[188] = 1000
p[189] = 1000
p[190] = 30
p[191] = 1000
p[192] = 1000
p[193] = 1000
p[194] = 1000
p[195] = 1000
p[196] = 1000
p[197] = 1000
p[198] = 3000
p[199] = 1000
p[200] = 1000
p[201] = 1000
p[202] = 100
p[203] = 100
p[204] = 100
p[205] = 100
p[206] = 1000
p[207] = 1000
p[208] = 1000
p[209] = 1000
p[210] = 1000
p[211] = 1000
p[212] = 1000
p[213] = 1000
p[214] = 1000
p[215] = 3000
p[216] = 100000
p[217] = 1000
p[218] = 1000
p[219] = 1000
p[220] = 1000
p[221] = 1000
p[222] = 1000
p[223] = 1000
p[224] = 1000
p[225] = 1000
p[226] = 1000
p[227] = 1000
p[228] = 0
p[229] = 1
p[230] = 0
p[231] = 0
p[232] = 1
p[233] = 1
p[234] = 1
p[235] = -0.7
p[236] = -0.7
p[237] = 0.7
p[238] = 0.7
p[239] = -0.2
p[240] = -0.2
p[241] = 1
p[242] = 1
p[243] = -0.7
p[244] = -0.7
p[245] = 0.05
p[246] = 1
p[247] = 0.5
p[248] = -0.7
p[249] = -0.5
p[250] = 2
p[251] = 0.7
p[252] = 0.7
p[253] = -2
p[254] = -0.2
p[255] = -0.2
p[256] = 1
p[257] = 0.7
p[258] = -0.2
p[259] = -0.21
p[260] = 0
p[261] = 1
p[262] = 0.7
p[263] = 1
p[264] = -0.2
p[265] = -0.2
p[266] = 0.05
p[267] = 0
p[268] = -0.01
p[269] = 0.7
p[270] = 1
p[271] = -0.5
p[272] = -1
p[273] = 0.04
p[274] = -0.01
p[275] = -0.01
p[276] = -0.01
p[277] = 1
p[278] = 1
p[279] = 0.5
p[280] = -1
p[281] = -1
p[282] = -0.5
p[283] = 0.1
p[284] = -0.1
p[285] = 0.5
p[286] = 1
p[287] = 1
p[288] = -1
p[289] = -0.3
p[290] = -1
p[291] = 1
p[292] = 0.7
p[293] = -0.4
p[294] = -0.7
p[295] = 1
p[296] = -0.7
p[297] = 0.7
p[298] = 2
p[299] = -0.2
p[300] = -2
p[301] = 1
p[302] = 0.7
p[303] = 0.7
p[304] = 0.2
p[305] = 0.7
p[306] = -0.5
p[307] = -0.21
p[308] = -0.2
p[309] = -0.2
p[310] = 0
p[311] = 0
p[312] = -0.01
p[313] = 0.7
p[314] = -0.2
p[315] = 2
p[316] = 2
p[317] = -2
p[318] = -2
p[319] = 2
p[320] = 2
p[321] = -2
p[322] = -2
p[323] = 1
p[324] = 1
p[325] = -1
p[326] = -1
p[327] = 1
p[328] = -0.5
p[329] = 0.7
p[330] = 0.7
p[331] = 2
p[332] = -0.2
p[333] = -2
p[334] = 1
p[335] = 1
p[336] = -0.7
p[337] = 1
p[338] = -0.7
p[339] = 1
p[340] = -0.7
p[341] = 1
p[342] = -0.7
p[343] = 2
p[344] = -2
p[345] = 0.7
p[346] = -0.2
p[347] = -0.2
for i in 1:length(par)
p[i] = par[i]
end
#endif /* FIXED */
#ifdef ASSIGNMENT
y = zeros(102)
y[1] = log(x[1])
y[2] = log(x[2])
y[3] = log(x[3])
y[4] = log(x[4])
y[5] = log(x[5])
y[6] = log(x[6])
y[7] = log(x[7])
y[8] = log(x[8])
y[9] = log(x[9])
y[10] = log(x[10])
y[11] = log(x[11])
y[12] = log(x[12])
y[13] = log(x[13])
y[14] = log(x[14])
y[15] = log(x[15])
y[16] = log(x[16])
y[17] = log(x[17])
y[18] = log(x[18])
y[19] = log(x[19])
y[20] = log(x[20])
y[21] = log(x[21])
y[22] = log(x[22])
y[23] = log(x[23])
y[24] = log(x[24])
y[25] = log(x[25])
y[26] = log(x[26])
y[27] = log(x[27])
y[28] = log(x[28])
y[29] = log(x[29])
y[30] = log(x[30])
y[31] = log(x[31])
y[32] = log(x[32])
y[33] = log(x[33])
y[34] = log(x[34])
y[35] = log(x[35])
y[36] = x[1]*p[186]
y[37] = x[2]*p[187]
y[38] = x[3]*p[188]
y[39] = x[4]*p[189]
y[40] = x[5]*p[190]
y[41] = x[6]*p[191]
y[42] = x[7]*p[192]
y[43] = x[8]*p[193]
y[44] = x[9]*p[194]
y[45] = x[10]*p[195]
y[46] = x[11]*p[196]
y[47] = x[12]*p[197]
y[48] = x[13]*p[198]
y[49] = x[14]*p[199]
y[50] = x[15]*p[200]
y[51] = x[16]*p[201]
y[52] = x[17]*p[202]
y[53] = x[18]*p[203]
y[54] = x[19]*p[204]
y[55] = x[20]*p[205]
y[56] = x[21]*p[206]
y[57] = x[22]*p[207]
y[58] = x[23]*p[208]
y[59] = x[24]*p[209]
y[60] = x[25]*p[210]
y[61] = x[26]*p[211]
y[62] = x[27]*p[212]
y[63] = x[28]*p[213]
y[64] = x[29]*p[214]
y[65] = x[30]*p[215]
y[66] = x[31]*p[216]
y[67] = x[32]*p[217]
y[68] = x[33]*p[218]
y[69] = x[34]*p[219]
if t>5 && t<10
y[70] = p[150]
else
y[70] = 0.00000000000000000
end
y[71] = 0.00000000000000000
y[72] = 0.00000000000000000
y[73] = 0.00000000000000000
y[74] = 0.00000000000000000
y[75] = 0.00000000000000000
y[76] = p[123]*p[159]*(1.00000000000000000+p[233]*y[8]+p[234]*y[9]+p[235]*y[6]+p[236]*y[7])
y[77] = p[124]*p[160]*(1.00000000000000000+p[237]*y[12]+p[238]*y[13]+p[239]*y[10]+p[240]*y[11])
y[78] = p[125]*p[161]*(1.00000000000000000+p[241]*y[15]+p[242]*y[6]+p[243]*y[8]+p[244]*y[14]+p[245]*y[16])
y[79] = p[126]*p[162]*(1.00000000000000000+p[246]*y[16]+p[247]*y[7]+p[248]*y[15]+p[249]*y[9])
y[80] = p[127]*p[163]*(1.00000000000000000+p[250]*y[9]+p[251]*y[19]+p[252]*y[20]+p[253]*y[7]+p[254]*y[17]+p[255]*y[18])
y[81] = p[128]*p[164]*(1.00000000000000000+p[256]*y[22]+p[257]*y[11]+p[258]*y[21]+p[259]*y[13]+p[260]*y[13])
y[82] = p[129]*p[165]*(1.00000000000000000+p[261]*y[15]+p[262]*y[7]+p[263]*y[21]+p[264]*y[8]+p[265]*y[9]+p[266]*y[32]+p[267]*y[9]+p[268]*y[30])
y[83] = p[130]*p[166]*(1.00000000000000000+p[269]*y[25]+p[270]*y[26]+p[271]*y[23]+p[272]*y[24]+p[273]*y[11]+p[274]*y[11]+p[275]*y[13]+p[276]*y[22])
y[84] = p[131]*p[167]*(1.00000000000000000+p[277]*y[28]+p[278]*y[29]+p[279]*y[23]+p[280]*y[12]+p[281]*y[27]+p[282]*y[25]+p[283]*y[27]+p[284]*y[10])
y[85] = p[132]*p[168]*(1.00000000000000000+p[285]*y[32]+p[286]*y[22]+p[287]*y[16]+p[288]*y[15]+p[289]*y[30]+p[290]*y[27])
y[86] = p[133]*p[169]*(1.00000000000000000+p[291]*y[21]+p[292]*y[23]+p[293]*y[2]+p[294]*y[25])
y[87] = p[134]*p[170]*(1.00000000000000000+p[295]*y[24]+p[296]*y[22])
y[88] = p[135]*p[171]*(1.00000000000000000+p[297]*y[20]+p[298]*y[17]+p[299]*y[18]+p[300]*y[19])
y[89] = p[136]*(p[172]/(p[209]/p[220]*(p[210]/p[221])*(p[197]/p[222])*(p[218]/p[223])*(p[219]/p[224])*(p[214]/p[225])*(p[195]/p[226])*(p[198]/p[227])/((1.00000000000000000+p[209]/p[220])*(1.00000000000000000+p[210]/p[221])*(1.00000000000000000+p[197]/p[222])*(1.00000000000000000+p[218]/p[223])*(1.00000000000000000+p[219]/p[224])*(1.00000000000000000+p[214]/p[225])*(1.00000000000000000+p[195]/p[226])*(1.00000000000000000+p[198]/p[227]))))*(x[24]*p[209]/p[220]*(x[25]*p[210]/p[221])*(x[12]*p[197]/p[222])*(x[33]*p[218]/p[223])*(x[34]*p[219]/p[224])*(x[29]*p[214]/p[225])*(x[10]*p[195]/p[226])*(x[13]*p[198]/p[227])/((1.00000000000000000+x[24]*p[209]/p[220])*(1.00000000000000000+x[25]*p[210]/p[221])*(1.00000000000000000+x[12]*p[197]/p[222])*(1.00000000000000000+x[33]*p[218]/p[223])*(1.00000000000000000+x[34]*p[219]/p[224])*(1.00000000000000000+x[29]*p[214]/p[225])*(1.00000000000000000+x[10]*p[195]/p[226])*(1.00000000000000000+x[13]*p[198]/p[227])))
y[90] = p[137]*p[173]*(1.00000000000000000+p[301]*y[8]+p[302]*y[7]+p[303]*y[19]+p[304]*y[31]+p[305]*y[32]+p[306]*y[16]+p[307]*y[9]+p[308]*y[17]+p[309]*y[30]+p[310]*y[31]+p[311]*y[9]+p[312]*y[15])
y[91] = p[138]*p[174]*(1.00000000000000000+p[313]*y[18]+p[314]*y[20])
y[92] = p[139]*p[175]*(1.00000000000000000+p[315]*y[11]+p[316]*y[30]+p[317]*y[32]+p[318]*y[13])
y[93] = p[140]*p[176]*(1.00000000000000000+p[319]*y[8]+p[320]*y[27]+p[321]*y[28]+p[322]*y[16])
y[94] = p[141]*p[177]*(1.00000000000000000+p[323]*y[14]+p[324]*y[12]+p[325]*y[29]+p[326]*y[6])
y[95] = p[142]*p[178]*(1.00000000000000000+p[327]*y[13]+p[328]*y[11])
y[96] = p[143]*p[179]*(1.00000000000000000+p[329]*y[32]+p[330]*y[31]+p[331]*y[18]+p[332]*y[30]+p[333]*y[20])
y[97] = p[144]*p[180]*(1.00000000000000000+p[334]*y[27]+p[335]*y[31]+p[336]*y[16])
y[98] = p[145]*p[181]*(1.00000000000000000+p[337]*y[1]+p[338]*y[26])
y[99] = p[146]*p[182]*(1.00000000000000000+p[339]*y[3]+p[340]*y[33])
y[100] = p[147]*p[183]*(1.00000000000000000+p[341]*y[4]+p[342]*y[34])
y[101] = p[148]*p[184]*(1.00000000000000000+p[343]*y[6]+p[344]*y[12])
y[102] = p[149]*p[185]*(1.00000000000000000+p[345]*y[18]+p[346]*y[31]+p[347]*y[20])
#endif /* ASSIGNMENT */
#ifdef ODEs
dx_dt[1] = (y[71]-y[98]*1.00000000000000000/p[151]-x[1]*p[186]/x[35]*y[70])/p[186]
dx_dt[2] = ((-y[72])+y[86]*1.00000000000000000/p[151]-x[2]*p[187]/x[35]*y[70])/p[187]
dx_dt[3] = (y[73]-y[99]*1.00000000000000000/p[151]-x[3]*p[188]/x[35]*y[70])/p[188]
dx_dt[4] = (y[74]-y[100]*1.00000000000000000/p[151]-x[4]*p[189]/x[35]*y[70])/p[189]
dx_dt[5] = ((-y[75])+y[89]*1.00000000000000000/p[151]-x[5]*p[190]/x[35]*y[70])/p[190]
dx_dt[6] = (y[76]-y[78]+y[94]*p[152]/p[153]-y[101]*p[152]/p[153])/p[191]
dx_dt[7] = (y[76]-y[79]+y[80]-2.00000000000000000*y[82]*p[152]/p[153]-y[90])/p[192]
dx_dt[8] = ((-y[76])+y[78]+y[82]*p[152]/p[153]-y[90]-y[93]*p[152]/p[153])/p[193]
dx_dt[9] = ((-y[76])+y[79]-y[80]+2.00000000000000000*y[82]*p[152]/p[153]+y[90])/p[194]
dx_dt[10] = (y[77]-120.00000000000000000*y[89])/p[195]
dx_dt[11] = (y[77]-y[81]+1260.00000000000000000*y[89]-y[92]+y[95])/p[196]
dx_dt[12] = ((-y[77])+y[84]-240.00000000000000000*y[89]-y[94]+y[101])/p[197]
dx_dt[13] = ((-y[77])+y[81]-1260.00000000000000000*y[89]+y[92]-y[95])/p[198]
dx_dt[14] = (y[78]-y[94]*p[152]/p[153])/p[199]
dx_dt[15] = ((-y[78])+y[79]-y[82]*p[152]/p[153]+y[85]*p[152]/p[153])/p[200]
dx_dt[16] = ((-y[79])-y[85]*p[152]/p[153]+y[90]+y[93]*p[152]/p[153]+y[97]*p[152]/p[153])/p[201]
dx_dt[17] = (2.00000000000000000*y[80]-2.00000000000000000*y[88]+2.00000000000000000*y[90])/p[202]
dx_dt[18] = (4.00000000000000000*y[80]+6.00000000000000000*y[88]-y[91]-3.00000000000000000*y[96]-y[102]*p[152]/p[153])/p[203]
dx_dt[19] = ((-2.00000000000000000)*y[80]+2.00000000000000000*y[88]-2.00000000000000000*y[90])/p[204]
dx_dt[20] = ((-4.00000000000000000)*y[80]-6.00000000000000000*y[88]+y[91]+3.00000000000000000*y[96]+y[102]*p[152]/p[153])/p[205]
dx_dt[21] = (y[81]-y[82]-y[86])/p[206]
dx_dt[22] = ((-y[81])-y[85]+y[87])/p[207]
dx_dt[23] = (y[83]-y[84]-y[86]+120.00000000000000000*y[89])/p[208]
dx_dt[24] = (y[83]-y[87]-120.00000000000000000*y[89])/p[209]
dx_dt[25] = ((-y[83])+y[84]+y[86]-120.00000000000000000*y[89])/p[210]
dx_dt[26] = ((-0.50000000000000000)*y[83]+y[98])/p[211]
dx_dt[27] = (y[84]+y[85]-y[93]-y[97])/p[212]
dx_dt[28] = ((-y[84])+120.00000000000000000*y[89]+y[93])/p[213]
dx_dt[29] = ((-y[84])-120.00000000000000000*y[89]+y[94])/p[214]
dx_dt[30] = (y[85]*p[152]/p[153]+y[90]-y[92]*p[152]/p[153]+y[96])/p[215]
dx_dt[31] = ((-y[90])-y[96]-y[97]*p[152]/p[153]+y[102]*p[152]/p[153])/p[216]
dx_dt[32] = ((-y[85])*p[152]/p[153]-y[90]+y[92]*p[152]/p[153]-y[96])/p[217]
dx_dt[33] = ((-120.00000000000000000)*y[89]+y[99])/p[218]
dx_dt[34] = ((-120.00000000000000000)*y[89]+y[100])/p[219]
dx_dt[35] = y[70]
end
tspan = (1.0, 20.0)
solvers_list = [lsoda(), CVODE_BDF(), Rosenbrock23(), TRBDF2(), ABDF2(), Rodas5(), Rodas4P(), Kvaerno5(), KenCarp4(), RadauIIA5()]
for solver in solvers_list
try
prob = ODEProblem(fun_b4, x0, tspan, p)
sol = solve(prob, solver)
plot_sol = plot(sol, title="$solver"[1:match(r"(\{)|(\()", "$solver").offset-1])
display(plot_sol)
catch e
@warn("\nException:\n$e\nFor solver:\n$solver\n")
end
end
#=>
# One plot for each chemical element
sol_u = reduce(hcat, sol.u)
for i in 1:size(sol_u)[1]
p = plot(sol.t, sol_u[i,:])
display(p)
end
<=#
And the errors:
┌ Warning:
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ Rosenbrock23{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), Val{:central})
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning:
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ TRBDF2{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},NLNewton{Nothing,Nothing},DataType,Float64,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), NLNewton{Nothing,Nothing}(nothing, nothing, 10), Val{:central}, true, :linear, 0.001)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning:
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ ABDF2{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},NLNewton{Nothing,Nothing},DataType,Nothing,Nothing,Float64,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), NLNewton{Nothing,Nothing}(nothing, nothing, 10), Val{:central}, nothing, nothing, true, :linear, 0.001)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning:
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ Rodas5{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), Val{:central})
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning:
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ Rodas4P{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), Val{:central})
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning:
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ Kvaerno5{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},NLNewton{Nothing,Nothing},DataType,Float64,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), NLNewton{Nothing,Nothing}(nothing, nothing, 10), Val{:central}, true, :linear, 0.001)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning:
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ KenCarp4{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},NLNewton{Nothing,Nothing},DataType,Float64,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), NLNewton{Nothing,Nothing}(nothing, nothing, 10), Val{:central}, true, :linear, 0.001)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
┌ Warning:
│ Exception:
│ MethodError(Float64, (Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(fun_b4),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64}}(1.6094379124341003,0.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0),), 0x0000000000006444)
│ For solver:
│ RadauIIA5{0,true,LinSolveFactorize{typeof(LinearAlgebra.lu!)},DataType,Float64,Nothing,:Predictive}(LinSolveFactorize{typeof(LinearAlgebra.lu!)}(LinearAlgebra.lu!, nothing), Val{:central}, true, :dense, 0.001, nothing, nothing, 10, 1)
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:567
Not sure if this is an actual issue or if I’m doing something wrong here.
As a side note: in case anyone knows a better way for plotting each chemical element, as mentioned in the last comment of my code, please let me know. 
– Edit:
I just ran this code again and this time lsoda got lost as well.
┌ Warning:
│ Exception:
│ DomainError(-0.14279345102360447, "log will only return a complex result if called with a complex argument. Try log(Complex(x)).")
│ For solver:
│ lsoda()
└ @ Main ~/git/ChemParamEst/replicate_error_b4.jl:561
So the error here is clearer for me, it’s having trouble with the log function.
But it seems like most of the time lsoda will get it right, with this error popping up just a couple of times after some 20+ runs. CVODE_BDF doesn’t seem to have a problem with log at all.