search

SciML / DifferentialEquations.jl

Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equations (SDEs), delay differential equations (DDEs), differential-algebraic equations (DAEs), and more in Julia.
https://docs.sciml.ai/DiffEqDocs/stable/
Other
2.85k stars 225 forks source link

Solve extremely slow on first run #303

Closed nealde closed 6 years ago

nealde commented 6 years ago

Hello,

I'm solving some stiff DAEs with a mass-matrix solver in the following code:

using DifferentialEquations
include("model.jl")
pars = [7.33234842e-08,6.40146313e-12,1.56954148e-13,1.65301662e-06
,1.28277079e-05,1.39115369e+00,1.70413971e+00,1.93944836e+00
,1.15245117e+03,3.46144727e+04,4.02951641e+04,3.24118063e-02
,2.70097144e-02,5.20209372e-01,3.74625146e-01,4.16417062e-01
,5.51529622e-09,6.95793778e-09,4.91284482e-05,4.88584483e-05
,1.28054144e-05,1.38146591e+02,9.15019798e+00,9.78830695e-01
,4.58723247e-01,4.32648301e-01,0.02]
y0 = yyy(pars)
@time P2D2c(zeros(n),y0,pars,0.) # 1 second
@time P2D2c(zeros(n),y0,pars,0.) # .000029 seconds
#function condition(u,t,integrator)
#    u[131]-3.0
#end
#affect!(integrator) = terminate!(integrator)
#cb = ContinuousCallback(condition,affect!)
tspan = (0.0,1600.)
prob = ODEProblem{true}(P2D2c,y0,tspan,pars,mass_matrix=MM)
@time sol = solve(prob, Rodas4P(), reltol=1e-6, force_dtmin=true) # 500 seconds??
@time sol = solve(prob, Rodas4P(), reltol=1e-6, force_dtmin=true) # 1-2 seconds

and for smaller systems (even a similar system with fewer equations, ~26), compilation takes just a few seconds, around 10 seconds with .01 seconds of execution time. However, for systems which are slightly larger, I get extremely long first-run times relative to subsequent execution times:

519.003612 seconds (664.19 M allocations: 18.791 GiB, 0.98% gc time) 2.769820 seconds (1.96 M allocations: 170.898 MiB, 8.22% gc time)

This doesn't seem like normal behavior. The calculations are all in-place, and the functions are in the attached file. This happens both in Linux and Windows with the following package status:

 - DiffEqBase                    3.6.2
 - DiffEqBiological              2.0.0
 - DiffEqCallbacks               1.0.0
 - DiffEqDevTools                1.0.0
 - DiffEqFinancial               1.0.0
 - DiffEqJump                    1.1.0
 - DiffEqMonteCarlo              0.10.0
 - DiffEqNoiseProcess            1.0.0
 - DiffEqPDEBase                 0.4.0
 - DiffEqParamEstim              1.0.0
 - DiffEqPhysics                 1.0.0
 - DiffEqSensitivity             1.0.1
 - DiffEqUncertainty             0.1.0

model.jl can be found below:

P2D2c = function(dy,y,pars,t)
    Dbulk = pars[1]
    Dsn = pars[2]
    Dsp = pars[3]
    F = 96487.000000
    R = 8.314300
    Rpn = pars[4]
    Rpp = pars[5]
    Tr = 298.150000
    brugn = pars[6]
    brugp = pars[7]
    brugs = pars[8]
    c0 = pars[9]
    ctn = pars[10]
    ctp = pars[11]
    efn = pars[12]
    efp = pars[13]
    en = pars[14]
    ep = pars[15]
    es = pars[16]
    kn = pars[17]
    kp = pars[18]
    ln1 = pars[19]
    lp = pars[20]
    ls = pars[21]
    sigmn = pars[22]
    sigmp = pars[23]
    socn = pars[24]
    socp = pars[25]
    t1 = pars[26]

    current = -1.0
    capacity=17.5
    rate=-2.
    iapp=capacity*rate
    up=[0.,0.]
    Kr = pars[27]

    ap = 3/Rpp*(1-ep-efp)

    an = 3/Rpn*(1-en-efn)

    sigmap = sigmp*(1-ep-efp)

    sigman = sigmn*(1-en-efn)

    D2pos = ep^brugp*Dbulk

    D2sep = es^brugs*Dbulk

    D2neg = en^brugn*Dbulk

    dy[1] = 1/ep*(64.*D2pos/lp^2*(y[60]-2.*y[1]+y[2])+2.*ap*(1.-1.*t1)*kp*(y[1]*c0)^.5*(-1.*y[110]*ctp+ctp)^.5*(y[110]*ctp)^.5*sinh(.5*F/R/Tr*(y[86]-1.*y[65]-1.*ocp_cathode(y[110],up)))/c0)
    dy[2] = 1/ep*(64.*D2pos/lp^2*(y[1]-2.*y[2]+y[3])+2.*ap*(1.-1.*t1)*kp*(y[2]*c0)^.5*(-1.*y[111]*ctp+ctp)^.5*(y[111]*ctp)^.5*sinh(.5*F/R/Tr*(y[87]-1.*y[66]-1.*ocp_cathode(y[111],up)))/c0)
    dy[3] = 1/ep*(64.*D2pos/lp^2*(y[2]-2.*y[3]+y[4])+2.*ap*(1.-1.*t1)*kp*(y[3]*c0)^.5*(-1.*y[112]*ctp+ctp)^.5*(y[112]*ctp)^.5*sinh(.5*F/R/Tr*(y[88]-1.*y[67]-1.*ocp_cathode(y[112],up)))/c0)
    dy[4] = 1/ep*(64.*D2pos/lp^2*(y[3]-2.*y[4]+y[5])+2.*ap*(1.-1.*t1)*kp*(y[4]*c0)^.5*(-1.*y[113]*ctp+ctp)^.5*(y[113]*ctp)^.5*sinh(.5*F/R/Tr*(y[89]-1.*y[68]-1.*ocp_cathode(y[113],up)))/c0)
    dy[5] = 1/ep*(64.*D2pos/lp^2*(y[4]-2.*y[5]+y[6])+2.*ap*(1.-1.*t1)*kp*(y[5]*c0)^.5*(-1.*y[114]*ctp+ctp)^.5*(y[114]*ctp)^.5*sinh(.5*F/R/Tr*(y[90]-1.*y[69]-1.*ocp_cathode(y[114],up)))/c0)
    dy[6] = 1/ep*(64.*D2pos/lp^2*(y[5]-2.*y[6]+y[7])+2.*ap*(1.-1.*t1)*kp*(y[6]*c0)^.5*(-1.*y[115]*ctp+ctp)^.5*(y[115]*ctp)^.5*sinh(.5*F/R/Tr*(y[91]-1.*y[70]-1.*ocp_cathode(y[115],up)))/c0)
    dy[7] = 1/ep*(64.*D2pos/lp^2*(y[6]-2.*y[7]+y[61])+2.*ap*(1.-1.*t1)*kp*(y[7]*c0)^.5*(-1.*y[116]*ctp+ctp)^.5*(y[116]*ctp)^.5*sinh(.5*F/R/Tr*(y[92]-1.*y[71]-1.*ocp_cathode(y[116],up)))/c0)
    dy[8] = 16./es*D2sep/ls^2*(y[61]-2.*y[8]+y[9])
    dy[9] = 16./es*D2sep/ls^2*(y[8]-2.*y[9]+y[10])
    dy[10] = 16./es*D2sep/ls^2*(y[9]-2.*y[10]+y[62])
    dy[11] = 1/en*(64.*D2neg/ln1^2*(y[62]-2.*y[11]+y[12])+2.*an*(1.-1.*t1)*kn*(y[11]*c0)^.5*(-1.*y[124]*ctn+ctn)^.5*(y[124]*ctn)^.5*sinh(.5*F/R/Tr*(y[95]-1.*y[77]+.5e-1-2.325*exp(-100.*y[124]^1.15)+.1721*tanh(20.000000*y[124]-21.000000)+.25e-2*tanh(39.347962*y[124]-37.241379)+.34e-1*tanh(89.411765*y[124]-6.0294118)+.2e-2*tanh(7.0422535*y[124]-1.3661972)+.155e-1*tanh(77.519380*y[124]-8.1395349)))/c0)
    dy[12] = 1/en*(64.*D2neg/ln1^2*(y[11]-2.*y[12]+y[13])+2.*an*(1.-1.*t1)*kn*(y[12]*c0)^.5*(-1.*y[125]*ctn+ctn)^.5*(y[125]*ctn)^.5*sinh(.5*F/R/Tr*(y[96]-1.*y[78]+.5e-1-2.325*exp(-100.*y[125]^1.15)+.1721*tanh(20.000000*y[125]-21.000000)+.25e-2*tanh(39.347962*y[125]-37.241379)+.34e-1*tanh(89.411765*y[125]-6.0294118)+.2e-2*tanh(7.0422535*y[125]-1.3661972)+.155e-1*tanh(77.519380*y[125]-8.1395349)))/c0)
    dy[13] = 1/en*(64.*D2neg/ln1^2*(y[12]-2.*y[13]+y[14])+2.*an*(1.-1.*t1)*kn*(y[13]*c0)^.5*(-1.*y[126]*ctn+ctn)^.5*(y[126]*ctn)^.5*sinh(.5*F/R/Tr*(y[97]-1.*y[79]+.5e-1-2.325*exp(-100.*y[126]^1.15)+.1721*tanh(20.000000*y[126]-21.000000)+.25e-2*tanh(39.347962*y[126]-37.241379)+.34e-1*tanh(89.411765*y[126]-6.0294118)+.2e-2*tanh(7.0422535*y[126]-1.3661972)+.155e-1*tanh(77.519380*y[126]-8.1395349)))/c0)
    dy[14] = 1/en*(64.*D2neg/ln1^2*(y[13]-2.*y[14]+y[15])+2.*an*(1.-1.*t1)*kn*(y[14]*c0)^.5*(-1.*y[127]*ctn+ctn)^.5*(y[127]*ctn)^.5*sinh(.5*F/R/Tr*(y[98]-1.*y[80]+.5e-1-2.325*exp(-100.*y[127]^1.15)+.1721*tanh(20.000000*y[127]-21.000000)+.25e-2*tanh(39.347962*y[127]-37.241379)+.34e-1*tanh(89.411765*y[127]-6.0294118)+.2e-2*tanh(7.0422535*y[127]-1.3661972)+.155e-1*tanh(77.519380*y[127]-8.1395349)))/c0)
    dy[15] = 1/en*(64.*D2neg/ln1^2*(y[14]-2.*y[15]+y[16])+2.*an*(1.-1.*t1)*kn*(y[15]*c0)^.5*(-1.*y[128]*ctn+ctn)^.5*(y[128]*ctn)^.5*sinh(.5*F/R/Tr*(y[99]-1.*y[81]+.5e-1-2.325*exp(-100.*y[128]^1.15)+.1721*tanh(20.000000*y[128]-21.000000)+.25e-2*tanh(39.347962*y[128]-37.241379)+.34e-1*tanh(89.411765*y[128]-6.0294118)+.2e-2*tanh(7.0422535*y[128]-1.3661972)+.155e-1*tanh(77.519380*y[128]-8.1395349)))/c0)
    dy[16] = 1/en*(64.*D2neg/ln1^2*(y[15]-2.*y[16]+y[17])+2.*an*(1.-1.*t1)*kn*(y[16]*c0)^.5*(-1.*y[129]*ctn+ctn)^.5*(y[129]*ctn)^.5*sinh(.5*F/R/Tr*(y[100]-1.*y[82]+.5e-1-2.325*exp(-100.*y[129]^1.15)+.1721*tanh(20.000000*y[129]-21.000000)+.25e-2*tanh(39.347962*y[129]-37.241379)+.34e-1*tanh(89.411765*y[129]-6.0294118)+.2e-2*tanh(7.0422535*y[129]-1.3661972)+.155e-1*tanh(77.519380*y[129]-8.1395349)))/c0)
    dy[17] = 1/en*(64.*D2neg/ln1^2*(y[16]-2.*y[17]+y[63])+2.*an*(1.-1.*t1)*kn*(y[17]*c0)^.5*(-1.*y[130]*ctn+ctn)^.5*(y[130]*ctn)^.5*sinh(.5*F/R/Tr*(y[101]-1.*y[83]+.5e-1-2.325*exp(-100.*y[130]^1.15)+.1721*tanh(20.000000*y[130]-21.000000)+.25e-2*tanh(39.347962*y[130]-37.241379)+.34e-1*tanh(89.411765*y[130]-6.0294118)+.2e-2*tanh(7.0422535*y[130]-1.3661972)+.155e-1*tanh(77.519380*y[130]-8.1395349)))/c0)
    dy[18] = 16.*(Dsp*(y[25]-1.*y[103])+Dsp*(y[103]-2.*y[18]+y[25]))/Rpp^2
    dy[19] = 16.*(Dsp*(y[26]-1.*y[104])+Dsp*(y[104]-2.*y[19]+y[26]))/Rpp^2
    dy[20] = 16.*(Dsp*(y[27]-1.*y[105])+Dsp*(y[105]-2.*y[20]+y[27]))/Rpp^2
    dy[21] = 16.*(Dsp*(y[28]-1.*y[106])+Dsp*(y[106]-2.*y[21]+y[28]))/Rpp^2
    dy[22] = 16.*(Dsp*(y[29]-1.*y[107])+Dsp*(y[107]-2.*y[22]+y[29]))/Rpp^2
    dy[23] = 16.*(Dsp*(y[30]-1.*y[108])+Dsp*(y[108]-2.*y[23]+y[30]))/Rpp^2
    dy[24] = 16.*(Dsp*(y[31]-1.*y[109])+Dsp*(y[109]-2.*y[24]+y[31]))/Rpp^2
    dy[25] = 4.*(2.*Dsp*(y[32]-1.*y[18])+4.*Dsp*(y[18]-2.*y[25]+y[32]))/Rpp^2
    dy[26] = 4.*(2.*Dsp*(y[33]-1.*y[19])+4.*Dsp*(y[19]-2.*y[26]+y[33]))/Rpp^2
    dy[27] = 4.*(2.*Dsp*(y[34]-1.*y[20])+4.*Dsp*(y[20]-2.*y[27]+y[34]))/Rpp^2
    dy[28] = 4.*(2.*Dsp*(y[35]-1.*y[21])+4.*Dsp*(y[21]-2.*y[28]+y[35]))/Rpp^2
    dy[29] = 4.*(2.*Dsp*(y[36]-1.*y[22])+4.*Dsp*(y[22]-2.*y[29]+y[36]))/Rpp^2
    dy[30] = 4.*(2.*Dsp*(y[37]-1.*y[23])+4.*Dsp*(y[23]-2.*y[30]+y[37]))/Rpp^2
    dy[31] = 4.*(2.*Dsp*(y[38]-1.*y[24])+4.*Dsp*(y[24]-2.*y[31]+y[38]))/Rpp^2
    dy[32] = 1.777777778*(3.*Dsp*(y[110]-1.*y[25])+9.*Dsp*(y[25]-2.*y[32]+y[110]))/Rpp^2
    dy[33] = 1.777777778*(3.*Dsp*(y[111]-1.*y[26])+9.*Dsp*(y[26]-2.*y[33]+y[111]))/Rpp^2
    dy[34] = 1.777777778*(3.*Dsp*(y[112]-1.*y[27])+9.*Dsp*(y[27]-2.*y[34]+y[112]))/Rpp^2
    dy[35] = 1.777777778*(3.*Dsp*(y[113]-1.*y[28])+9.*Dsp*(y[28]-2.*y[35]+y[113]))/Rpp^2
    dy[36] = 1.777777778*(3.*Dsp*(y[114]-1.*y[29])+9.*Dsp*(y[29]-2.*y[36]+y[114]))/Rpp^2
    dy[37] = 1.777777778*(3.*Dsp*(y[115]-1.*y[30])+9.*Dsp*(y[30]-2.*y[37]+y[115]))/Rpp^2
    dy[38] = 1.777777778*(3.*Dsp*(y[116]-1.*y[31])+9.*Dsp*(y[31]-2.*y[38]+y[116]))/Rpp^2
    dy[39] = 16.*Dsn*(2.*y[46]-2.*y[39])/Rpn^2
    dy[40] = 16.*Dsn*(2.*y[47]-2.*y[40])/Rpn^2
    dy[41] = 16.*Dsn*(2.*y[48]-2.*y[41])/Rpn^2
    dy[42] = 16.*Dsn*(2.*y[49]-2.*y[42])/Rpn^2
    dy[43] = 16.*Dsn*(2.*y[50]-2.*y[43])/Rpn^2
    dy[44] = 16.*Dsn*(2.*y[51]-2.*y[44])/Rpn^2
    dy[45] = 16.*Dsn*(2.*y[52]-2.*y[45])/Rpn^2
    dy[46] = 4.*Dsn*(6.*y[53]+2.*y[39]-8.*y[46])/Rpn^2
    dy[47] = 4.*Dsn*(6.*y[54]+2.*y[40]-8.*y[47])/Rpn^2
    dy[48] = 4.*Dsn*(6.*y[55]+2.*y[41]-8.*y[48])/Rpn^2
    dy[49] = 4.*Dsn*(6.*y[56]+2.*y[42]-8.*y[49])/Rpn^2
    dy[50] = 4.*Dsn*(6.*y[57]+2.*y[43]-8.*y[50])/Rpn^2
    dy[51] = 4.*Dsn*(6.*y[58]+2.*y[44]-8.*y[51])/Rpn^2
    dy[52] = 4.*Dsn*(6.*y[59]+2.*y[45]-8.*y[52])/Rpn^2
    dy[53] = 1.777777778*Dsn*(12.*y[124]+6.*y[46]-18.*y[53])/Rpn^2
    dy[54] = 1.777777778*Dsn*(12.*y[125]+6.*y[47]-18.*y[54])/Rpn^2
    dy[55] = 1.777777778*Dsn*(12.*y[126]+6.*y[48]-18.*y[55])/Rpn^2
    dy[56] = 1.777777778*Dsn*(12.*y[127]+6.*y[49]-18.*y[56])/Rpn^2
    dy[57] = 1.777777778*Dsn*(12.*y[128]+6.*y[50]-18.*y[57])/Rpn^2
    dy[58] = 1.777777778*Dsn*(12.*y[129]+6.*y[51]-18.*y[58])/Rpn^2
    dy[59] = 1.777777778*Dsn*(12.*y[130]+6.*y[52]-18.*y[59])/Rpn^2
    dy[60] = 4.*D2pos*(-1.*y[2]-3.*y[60]+4.*y[1])/lp
    dy[61] = 4.*D2pos*(y[6]+3.*y[61]-4.*y[7])/lp-2.*D2sep*(-1.*y[9]-3.*y[61]+4.*y[8])/ls
    dy[62] = 2.*D2sep*(y[9]+3.*y[62]-4.*y[10])/ls-4.*D2neg*(-1.*y[12]-3.*y[62]+4.*y[11])/ln1
    dy[63] = 4.*D2neg*(y[16]+3.*y[63]-4.*y[17])/ln1
    dy[64] = 4.*(-1.*y[66]-3.*y[64]+4.*y[65])/lp
    dy[65] = -4.*sigmap/lp*(y[87]-1.*y[85])-4.*ep^brugp*(.41253e-1+.5007e-3*y[1]-.47212e-6*y[1]^2+.15904e-9*y[1]^3-.16018e-13*y[1]^4)/lp*(y[66]-1.*y[64])+8.*ep^brugp*(.41253e-1+.5007e-3*y[1]-.47212e-6*y[1]^2+.15904e-9*y[1]^3-.16018e-13*y[1]^4)*R*Tr/F*(1.-1.*t1)/lp*(y[2]-1.*y[60])/y[1]-iapp
    dy[66] = -4.*sigmap/lp*(y[88]-1.*y[86])-4.*ep^brugp*(.41253e-1+.5007e-3*y[2]-.47212e-6*y[2]^2+.15904e-9*y[2]^3-.16018e-13*y[2]^4)/lp*(y[67]-1.*y[65])+8.*ep^brugp*(.41253e-1+.5007e-3*y[2]-.47212e-6*y[2]^2+.15904e-9*y[2]^3-.16018e-13*y[2]^4)*R*Tr/F*(1.-1.*t1)/lp*(y[3]-1.*y[1])/y[2]-iapp
    dy[67] = -4.*sigmap/lp*(y[89]-1.*y[87])-4.*ep^brugp*(.41253e-1+.5007e-3*y[3]-.47212e-6*y[3]^2+.15904e-9*y[3]^3-.16018e-13*y[3]^4)/lp*(y[68]-1.*y[66])+8.*ep^brugp*(.41253e-1+.5007e-3*y[3]-.47212e-6*y[3]^2+.15904e-9*y[3]^3-.16018e-13*y[3]^4)*R*Tr/F*(1.-1.*t1)/lp*(y[4]-1.*y[2])/y[3]-iapp
    dy[68] = -4.*sigmap/lp*(y[90]-1.*y[88])-4.*ep^brugp*(.41253e-1+.5007e-3*y[4]-.47212e-6*y[4]^2+.15904e-9*y[4]^3-.16018e-13*y[4]^4)/lp*(y[69]-1.*y[67])+8.*ep^brugp*(.41253e-1+.5007e-3*y[4]-.47212e-6*y[4]^2+.15904e-9*y[4]^3-.16018e-13*y[4]^4)*R*Tr/F*(1.-1.*t1)/lp*(y[5]-1.*y[3])/y[4]-iapp
    dy[69] = -4.*sigmap/lp*(y[91]-1.*y[89])-4.*ep^brugp*(.41253e-1+.5007e-3*y[5]-.47212e-6*y[5]^2+.15904e-9*y[5]^3-.16018e-13*y[5]^4)/lp*(y[70]-1.*y[68])+8.*ep^brugp*(.41253e-1+.5007e-3*y[5]-.47212e-6*y[5]^2+.15904e-9*y[5]^3-.16018e-13*y[5]^4)*R*Tr/F*(1.-1.*t1)/lp*(y[6]-1.*y[4])/y[5]-iapp
    dy[70] = -4.*sigmap/lp*(y[92]-1.*y[90])-4.*ep^brugp*(.41253e-1+.5007e-3*y[6]-.47212e-6*y[6]^2+.15904e-9*y[6]^3-.16018e-13*y[6]^4)/lp*(y[71]-1.*y[69])+8.*ep^brugp*(.41253e-1+.5007e-3*y[6]-.47212e-6*y[6]^2+.15904e-9*y[6]^3-.16018e-13*y[6]^4)*R*Tr/F*(1.-1.*t1)/lp*(y[7]-1.*y[5])/y[6]-iapp
    dy[71] = -4.*sigmap/lp*(y[93]-1.*y[91])-4.*ep^brugp*(.41253e-1+.5007e-3*y[7]-.47212e-6*y[7]^2+.15904e-9*y[7]^3-.16018e-13*y[7]^4)/lp*(y[72]-1.*y[70])+8.*ep^brugp*(.41253e-1+.5007e-3*y[7]-.47212e-6*y[7]^2+.15904e-9*y[7]^3-.16018e-13*y[7]^4)*R*Tr/F*(1.-1.*t1)/lp*(y[61]-1.*y[6])/y[7]-iapp
    dy[72] = 4.*ep^brugp*(.41253e-1+.5007e-3*y[61]-.47212e-6*y[61]^2+.15904e-9*y[61]^3-.16018e-13*y[61]^4)*(y[70]+3.*y[72]-4.*y[71])/lp-2.*es^brugs*(.41253e-1+.5007e-3*y[61]-.47212e-6*y[61]^2+.15904e-9*y[61]^3-.16018e-13*y[61]^4)*(-1.*y[74]-3.*y[72]+4.*y[73])/ls
    dy[73] = -2.*es^brugs*(.41253e-1+.5007e-3*y[8]-.47212e-6*y[8]^2+.15904e-9*y[8]^3-.16018e-13*y[8]^4)/ls*(y[74]-1.*y[72])+4.*es^brugs*(.41253e-1+.5007e-3*y[8]-.47212e-6*y[8]^2+.15904e-9*y[8]^3-.16018e-13*y[8]^4)*R*Tr/F*(1.-1.*t1)/ls*(y[9]-1.*y[61])/y[8]-iapp
    dy[74] = -2.*es^brugs*(.41253e-1+.5007e-3*y[9]-.47212e-6*y[9]^2+.15904e-9*y[9]^3-.16018e-13*y[9]^4)/ls*(y[75]-1.*y[73])+4.*es^brugs*(.41253e-1+.5007e-3*y[9]-.47212e-6*y[9]^2+.15904e-9*y[9]^3-.16018e-13*y[9]^4)*R*Tr/F*(1.-1.*t1)/ls*(y[10]-1.*y[8])/y[9]-iapp
    dy[75] = -2.*es^brugs*(.41253e-1+.5007e-3*y[10]-.47212e-6*y[10]^2+.15904e-9*y[10]^3-.16018e-13*y[10]^4)/ls*(y[76]-1.*y[74])+4.*es^brugs*(.41253e-1+.5007e-3*y[10]-.47212e-6*y[10]^2+.15904e-9*y[10]^3-.16018e-13*y[10]^4)*R*Tr/F*(1.-1.*t1)/ls*(y[62]-1.*y[9])/y[10]-iapp
    dy[76] = 2.*es^brugs*(.41253e-1+.5007e-3*y[62]-.47212e-6*y[62]^2+.15904e-9*y[62]^3-.16018e-13*y[62]^4)*(y[74]+3.*y[76]-4.*y[75])/ls-4.*en^brugn*(.41253e-1+.5007e-3*y[62]-.47212e-6*y[62]^2+.15904e-9*y[62]^3-.16018e-13*y[62]^4)*(-1.*y[78]-3.*y[76]+4.*y[77])/ln1
    dy[77] = -4.*sigman/ln1*(y[96]-1.*y[94])-4.*en^brugn*(.41253e-1+.5007e-3*y[11]-.47212e-6*y[11]^2+.15904e-9*y[11]^3-.16018e-13*y[11]^4)/ln1*(y[78]-1.*y[76])+8.*en^brugn*(.41253e-1+.5007e-3*y[11]-.47212e-6*y[11]^2+.15904e-9*y[11]^3-.16018e-13*y[11]^4)*R*Tr/F*(1.-1.*t1)/ln1*(y[12]-1.*y[62])/y[11]-iapp
    dy[78] = -4.*sigman/ln1*(y[97]-1.*y[95])-4.*en^brugn*(.41253e-1+.5007e-3*y[12]-.47212e-6*y[12]^2+.15904e-9*y[12]^3-.16018e-13*y[12]^4)/ln1*(y[79]-1.*y[77])+8.*en^brugn*(.41253e-1+.5007e-3*y[12]-.47212e-6*y[12]^2+.15904e-9*y[12]^3-.16018e-13*y[12]^4)*R*Tr/F*(1.-1.*t1)/ln1*(y[13]-1.*y[11])/y[12]-iapp
    dy[79] = -4.*sigman/ln1*(y[98]-1.*y[96])-4.*en^brugn*(.41253e-1+.5007e-3*y[13]-.47212e-6*y[13]^2+.15904e-9*y[13]^3-.16018e-13*y[13]^4)/ln1*(y[80]-1.*y[78])+8.*en^brugn*(.41253e-1+.5007e-3*y[13]-.47212e-6*y[13]^2+.15904e-9*y[13]^3-.16018e-13*y[13]^4)*R*Tr/F*(1.-1.*t1)/ln1*(y[14]-1.*y[12])/y[13]-iapp
    dy[80] = -4.*sigman/ln1*(y[99]-1.*y[97])-4.*en^brugn*(.41253e-1+.5007e-3*y[14]-.47212e-6*y[14]^2+.15904e-9*y[14]^3-.16018e-13*y[14]^4)/ln1*(y[81]-1.*y[79])+8.*en^brugn*(.41253e-1+.5007e-3*y[14]-.47212e-6*y[14]^2+.15904e-9*y[14]^3-.16018e-13*y[14]^4)*R*Tr/F*(1.-1.*t1)/ln1*(y[15]-1.*y[13])/y[14]-iapp
    dy[81] = -4.*sigman/ln1*(y[100]-1.*y[98])-4.*en^brugn*(.41253e-1+.5007e-3*y[15]-.47212e-6*y[15]^2+.15904e-9*y[15]^3-.16018e-13*y[15]^4)/ln1*(y[82]-1.*y[80])+8.*en^brugn*(.41253e-1+.5007e-3*y[15]-.47212e-6*y[15]^2+.15904e-9*y[15]^3-.16018e-13*y[15]^4)*R*Tr/F*(1.-1.*t1)/ln1*(y[16]-1.*y[14])/y[15]-iapp
    dy[82] = -4.*sigman/ln1*(y[101]-1.*y[99])-4.*en^brugn*(.41253e-1+.5007e-3*y[16]-.47212e-6*y[16]^2+.15904e-9*y[16]^3-.16018e-13*y[16]^4)/ln1*(y[83]-1.*y[81])+8.*en^brugn*(.41253e-1+.5007e-3*y[16]-.47212e-6*y[16]^2+.15904e-9*y[16]^3-.16018e-13*y[16]^4)*R*Tr/F*(1.-1.*t1)/ln1*(y[17]-1.*y[15])/y[16]-iapp
    dy[83] = -4.*sigman/ln1*(y[102]-1.*y[100])-4.*en^brugn*(.41253e-1+.5007e-3*y[17]-.47212e-6*y[17]^2+.15904e-9*y[17]^3-.16018e-13*y[17]^4)/ln1*(y[84]-1.*y[82])+8.*en^brugn*(.41253e-1+.5007e-3*y[17]-.47212e-6*y[17]^2+.15904e-9*y[17]^3-.16018e-13*y[17]^4)*R*Tr/F*(1.-1.*t1)/ln1*(y[63]-1.*y[16])/y[17]-iapp
    dy[84] = y[84]
    dy[85] = 4.*(-1.*y[87]-3.*y[85]+4.*y[86])/lp+1.*iapp/sigmap
    dy[86] = 64.*sigmap/lp^2*(y[85]-2.*y[86]+y[87])-2.*ap*F*kp*(y[1]*c0)^.5*(-1.*y[110]*ctp+ctp)^.5*(y[110]*ctp)^.5*sinh(.5*F/R/Tr*(y[86]-1.*y[65]-1.*ocp_cathode(y[110],up)))
    dy[87] = 64.*sigmap/lp^2*(y[86]-2.*y[87]+y[88])-2.*ap*F*kp*(y[2]*c0)^.5*(-1.*y[111]*ctp+ctp)^.5*(y[111]*ctp)^.5*sinh(.5*F/R/Tr*(y[87]-1.*y[66]-1.*ocp_cathode(y[111],up)))
    dy[88] = 64.*sigmap/lp^2*(y[87]-2.*y[88]+y[89])-2.*ap*F*kp*(y[3]*c0)^.5*(-1.*y[112]*ctp+ctp)^.5*(y[112]*ctp)^.5*sinh(.5*F/R/Tr*(y[88]-1.*y[67]-1.*ocp_cathode(y[112],up)))
    dy[89] = 64.*sigmap/lp^2*(y[88]-2.*y[89]+y[90])-2.*ap*F*kp*(y[4]*c0)^.5*(-1.*y[113]*ctp+ctp)^.5*(y[113]*ctp)^.5*sinh(.5*F/R/Tr*(y[89]-1.*y[68]-1.*ocp_cathode(y[113],up)))
    dy[90] = 64.*sigmap/lp^2*(y[89]-2.*y[90]+y[91])-2.*ap*F*kp*(y[5]*c0)^.5*(-1.*y[114]*ctp+ctp)^.5*(y[114]*ctp)^.5*sinh(.5*F/R/Tr*(y[90]-1.*y[69]-1.*ocp_cathode(y[114],up)))
    dy[91] = 64.*sigmap/lp^2*(y[90]-2.*y[91]+y[92])-2.*ap*F*kp*(y[6]*c0)^.5*(-1.*y[115]*ctp+ctp)^.5*(y[115]*ctp)^.5*sinh(.5*F/R/Tr*(y[91]-1.*y[70]-1.*ocp_cathode(y[115],up)))
    dy[92] = 64.*sigmap/lp^2*(y[91]-2.*y[92]+y[93])-2.*ap*F*kp*(y[7]*c0)^.5*(-1.*y[116]*ctp+ctp)^.5*(y[116]*ctp)^.5*sinh(.5*F/R/Tr*(y[92]-1.*y[71]-1.*ocp_cathode(y[116],up)))
    dy[93] = 4.*(y[91]+3.*y[93]-4.*y[92])/lp
    dy[94] = 4.*(-1.*y[96]-3.*y[94]+4.*y[95])/ln1
    dy[95] = 64.*sigman/ln1^2*(y[94]-2.*y[95]+y[96])-2.*an*F*kn*(y[11]*c0)^.5*(-1.*y[124]*ctn+ctn)^.5*(y[124]*ctn)^.5*sinh(.5*F/R/Tr*(y[95]-1.*y[77]+.5e-1-2.325*exp(-100.*y[124]^1.15)+.1721*tanh(20.000000*y[124]-21.000000)+.25e-2*tanh(39.347962*y[124]-37.241379)+.34e-1*tanh(89.411765*y[124]-6.0294118)+.2e-2*tanh(7.0422535*y[124]-1.3661972)+.155e-1*tanh(77.519380*y[124]-8.1395349)))
    dy[96] = 64.*sigman/ln1^2*(y[95]-2.*y[96]+y[97])-2.*an*F*kn*(y[12]*c0)^.5*(-1.*y[125]*ctn+ctn)^.5*(y[125]*ctn)^.5*sinh(.5*F/R/Tr*(y[96]-1.*y[78]+.5e-1-2.325*exp(-100.*y[125]^1.15)+.1721*tanh(20.000000*y[125]-21.000000)+.25e-2*tanh(39.347962*y[125]-37.241379)+.34e-1*tanh(89.411765*y[125]-6.0294118)+.2e-2*tanh(7.0422535*y[125]-1.3661972)+.155e-1*tanh(77.519380*y[125]-8.1395349)))
    dy[97] = 64.*sigman/ln1^2*(y[96]-2.*y[97]+y[98])-2.*an*F*kn*(y[13]*c0)^.5*(-1.*y[126]*ctn+ctn)^.5*(y[126]*ctn)^.5*sinh(.5*F/R/Tr*(y[97]-1.*y[79]+.5e-1-2.325*exp(-100.*y[126]^1.15)+.1721*tanh(20.000000*y[126]-21.000000)+.25e-2*tanh(39.347962*y[126]-37.241379)+.34e-1*tanh(89.411765*y[126]-6.0294118)+.2e-2*tanh(7.0422535*y[126]-1.3661972)+.155e-1*tanh(77.519380*y[126]-8.1395349)))
    dy[98] = 64.*sigman/ln1^2*(y[97]-2.*y[98]+y[99])-2.*an*F*kn*(y[14]*c0)^.5*(-1.*y[127]*ctn+ctn)^.5*(y[127]*ctn)^.5*sinh(.5*F/R/Tr*(y[98]-1.*y[80]+.5e-1-2.325*exp(-100.*y[127]^1.15)+.1721*tanh(20.000000*y[127]-21.000000)+.25e-2*tanh(39.347962*y[127]-37.241379)+.34e-1*tanh(89.411765*y[127]-6.0294118)+.2e-2*tanh(7.0422535*y[127]-1.3661972)+.155e-1*tanh(77.519380*y[127]-8.1395349)))
    dy[99] = 64.*sigman/ln1^2*(y[98]-2.*y[99]+y[100])-2.*an*F*kn*(y[15]*c0)^.5*(-1.*y[128]*ctn+ctn)^.5*(y[128]*ctn)^.5*sinh(.5*F/R/Tr*(y[99]-1.*y[81]+.5e-1-2.325*exp(-100.*y[128]^1.15)+.1721*tanh(20.000000*y[128]-21.000000)+.25e-2*tanh(39.347962*y[128]-37.241379)+.34e-1*tanh(89.411765*y[128]-6.0294118)+.2e-2*tanh(7.0422535*y[128]-1.3661972)+.155e-1*tanh(77.519380*y[128]-8.1395349)))
    dy[100] = 64.*sigman/ln1^2*(y[99]-2.*y[100]+y[101])-2.*an*F*kn*(y[16]*c0)^.5*(-1.*y[129]*ctn+ctn)^.5*(y[129]*ctn)^.5*sinh(.5*F/R/Tr*(y[100]-1.*y[82]+.5e-1-2.325*exp(-100.*y[129]^1.15)+.1721*tanh(20.000000*y[129]-21.000000)+.25e-2*tanh(39.347962*y[129]-37.241379)+.34e-1*tanh(89.411765*y[129]-6.0294118)+.2e-2*tanh(7.0422535*y[129]-1.3661972)+.155e-1*tanh(77.519380*y[129]-8.1395349)))
    dy[101] = 64.*sigman/ln1^2*(y[100]-2.*y[101]+y[102])-2.*an*F*kn*(y[17]*c0)^.5*(-1.*y[130]*ctn+ctn)^.5*(y[130]*ctn)^.5*sinh(.5*F/R/Tr*(y[101]-1.*y[83]+.5e-1-2.325*exp(-100.*y[130]^1.15)+.1721*tanh(20.000000*y[130]-21.000000)+.25e-2*tanh(39.347962*y[130]-37.241379)+.34e-1*tanh(89.411765*y[130]-6.0294118)+.2e-2*tanh(7.0422535*y[130]-1.3661972)+.155e-1*tanh(77.519380*y[130]-8.1395349)))
    dy[102] = 4.*(y[100]+3.*y[102]-4.*y[101])/ln1+1.*iapp/sigman
    dy[103] = -2.*(-1.*y[25]-3.*y[103]+4.*y[18])/Rpp
    dy[104] = -2.*(-1.*y[26]-3.*y[104]+4.*y[19])/Rpp
    dy[105] = -2.*(-1.*y[27]-3.*y[105]+4.*y[20])/Rpp
    dy[106] = -2.*(-1.*y[28]-3.*y[106]+4.*y[21])/Rpp
    dy[107] = -2.*(-1.*y[29]-3.*y[107]+4.*y[22])/Rpp
    dy[108] = -2.*(-1.*y[30]-3.*y[108]+4.*y[23])/Rpp
    dy[109] = -2.*(-1.*y[31]-3.*y[109]+4.*y[24])/Rpp
    dy[110] = 2.*Dsp*(y[25]+3.*y[110]-4.*y[32])/Rpp+2.*kp*(y[1]*c0)^.5*(-1.*y[110]*ctp+ctp)^.5*(y[110]*ctp)^.5*sinh(.5*F/R/Tr*(y[86]-1.*y[65]-1.*ocp_cathode(y[110],up)))/ctp
    dy[111] = 2.*Dsp*(y[26]+3.*y[111]-4.*y[33])/Rpp+2.*kp*(y[2]*c0)^.5*(-1.*y[111]*ctp+ctp)^.5*(y[111]*ctp)^.5*sinh(.5*F/R/Tr*(y[87]-1.*y[66]-1.*ocp_cathode(y[111],up)))/ctp
    dy[112] = 2.*Dsp*(y[27]+3.*y[112]-4.*y[34])/Rpp+2.*kp*(y[3]*c0)^.5*(-1.*y[112]*ctp+ctp)^.5*(y[112]*ctp)^.5*sinh(.5*F/R/Tr*(y[88]-1.*y[67]-1.*ocp_cathode(y[112],up)))/ctp
    dy[113] = 2.*Dsp*(y[28]+3.*y[113]-4.*y[35])/Rpp+2.*kp*(y[4]*c0)^.5*(-1.*y[113]*ctp+ctp)^.5*(y[113]*ctp)^.5*sinh(.5*F/R/Tr*(y[89]-1.*y[68]-1.*ocp_cathode(y[113],up)))/ctp
    dy[114] = 2.*Dsp*(y[29]+3.*y[114]-4.*y[36])/Rpp+2.*kp*(y[5]*c0)^.5*(-1.*y[114]*ctp+ctp)^.5*(y[114]*ctp)^.5*sinh(.5*F/R/Tr*(y[90]-1.*y[69]-1.*ocp_cathode(y[114],up)))/ctp
    dy[115] = 2.*Dsp*(y[30]+3.*y[115]-4.*y[37])/Rpp+2.*kp*(y[6]*c0)^.5*(-1.*y[115]*ctp+ctp)^.5*(y[115]*ctp)^.5*sinh(.5*F/R/Tr*(y[91]-1.*y[70]-1.*ocp_cathode(y[115],up)))/ctp
    dy[116] = 2.*Dsp*(y[31]+3.*y[116]-4.*y[38])/Rpp+2.*kp*(y[7]*c0)^.5*(-1.*y[116]*ctp+ctp)^.5*(y[116]*ctp)^.5*sinh(.5*F/R/Tr*(y[92]-1.*y[71]-1.*ocp_cathode(y[116],up)))/ctp
    dy[117] = -2.*(-1.*y[46]-3.*y[117]+4.*y[39])/Rpn
    dy[118] = -2.*(-1.*y[47]-3.*y[118]+4.*y[40])/Rpn
    dy[119] = -2.*(-1.*y[48]-3.*y[119]+4.*y[41])/Rpn
    dy[120] = -2.*(-1.*y[49]-3.*y[120]+4.*y[42])/Rpn
    dy[121] = -2.*(-1.*y[50]-3.*y[121]+4.*y[43])/Rpn
    dy[122] = -2.*(-1.*y[51]-3.*y[122]+4.*y[44])/Rpn
    dy[123] = -2.*(-1.*y[52]-3.*y[123]+4.*y[45])/Rpn
    dy[124] = 2.*Dsn*(y[46]+3.*y[124]-4.*y[53])/Rpn+2.*kn*(y[11]*c0)^.5*(-1.*y[124]*ctn+ctn)^.5*(y[124]*ctn)^.5*sinh(.5*F/R/Tr*(y[95]-1.*y[77]+.5e-1-2.325*exp(-100.*y[124]^1.15)+.1721*tanh(20.000000*y[124]-21.000000)+.25e-2*tanh(39.347962*y[124]-37.241379)+.34e-1*tanh(89.411765*y[124]-6.0294118)+.2e-2*tanh(7.0422535*y[124]-1.3661972)+.155e-1*tanh(77.519380*y[124]-8.1395349)))/ctn
    dy[125] = 2.*Dsn*(y[47]+3.*y[125]-4.*y[54])/Rpn+2.*kn*(y[12]*c0)^.5*(-1.*y[125]*ctn+ctn)^.5*(y[125]*ctn)^.5*sinh(.5*F/R/Tr*(y[96]-1.*y[78]+.5e-1-2.325*exp(-100.*y[125]^1.15)+.1721*tanh(20.000000*y[125]-21.000000)+.25e-2*tanh(39.347962*y[125]-37.241379)+.34e-1*tanh(89.411765*y[125]-6.0294118)+.2e-2*tanh(7.0422535*y[125]-1.3661972)+.155e-1*tanh(77.519380*y[125]-8.1395349)))/ctn
    dy[126] = 2.*Dsn*(y[48]+3.*y[126]-4.*y[55])/Rpn+2.*kn*(y[13]*c0)^.5*(-1.*y[126]*ctn+ctn)^.5*(y[126]*ctn)^.5*sinh(.5*F/R/Tr*(y[97]-1.*y[79]+.5e-1-2.325*exp(-100.*y[126]^1.15)+.1721*tanh(20.000000*y[126]-21.000000)+.25e-2*tanh(39.347962*y[126]-37.241379)+.34e-1*tanh(89.411765*y[126]-6.0294118)+.2e-2*tanh(7.0422535*y[126]-1.3661972)+.155e-1*tanh(77.519380*y[126]-8.1395349)))/ctn
    dy[127] = 2.*Dsn*(y[49]+3.*y[127]-4.*y[56])/Rpn+2.*kn*(y[14]*c0)^.5*(-1.*y[127]*ctn+ctn)^.5*(y[127]*ctn)^.5*sinh(.5*F/R/Tr*(y[98]-1.*y[80]+.5e-1-2.325*exp(-100.*y[127]^1.15)+.1721*tanh(20.000000*y[127]-21.000000)+.25e-2*tanh(39.347962*y[127]-37.241379)+.34e-1*tanh(89.411765*y[127]-6.0294118)+.2e-2*tanh(7.0422535*y[127]-1.3661972)+.155e-1*tanh(77.519380*y[127]-8.1395349)))/ctn
    dy[128] = 2.*Dsn*(y[50]+3.*y[128]-4.*y[57])/Rpn+2.*kn*(y[15]*c0)^.5*(-1.*y[128]*ctn+ctn)^.5*(y[128]*ctn)^.5*sinh(.5*F/R/Tr*(y[99]-1.*y[81]+.5e-1-2.325*exp(-100.*y[128]^1.15)+.1721*tanh(20.000000*y[128]-21.000000)+.25e-2*tanh(39.347962*y[128]-37.241379)+.34e-1*tanh(89.411765*y[128]-6.0294118)+.2e-2*tanh(7.0422535*y[128]-1.3661972)+.155e-1*tanh(77.519380*y[128]-8.1395349)))/ctn
    dy[129] = 2.*Dsn*(y[51]+3.*y[129]-4.*y[58])/Rpn+2.*kn*(y[16]*c0)^.5*(-1.*y[129]*ctn+ctn)^.5*(y[129]*ctn)^.5*sinh(.5*F/R/Tr*(y[100]-1.*y[82]+.5e-1-2.325*exp(-100.*y[129]^1.15)+.1721*tanh(20.000000*y[129]-21.000000)+.25e-2*tanh(39.347962*y[129]-37.241379)+.34e-1*tanh(89.411765*y[129]-6.0294118)+.2e-2*tanh(7.0422535*y[129]-1.3661972)+.155e-1*tanh(77.519380*y[129]-8.1395349)))/ctn
    dy[130] = 2.*Dsn*(y[52]+3.*y[130]-4.*y[59])/Rpn+2.*kn*(y[17]*c0)^.5*(-1.*y[130]*ctn+ctn)^.5*(y[130]*ctn)^.5*sinh(.5*F/R/Tr*(y[101]-1.*y[83]+.5e-1-2.325*exp(-100.*y[130]^1.15)+.1721*tanh(20.000000*y[130]-21.000000)+.25e-2*tanh(39.347962*y[130]-37.241379)+.34e-1*tanh(89.411765*y[130]-6.0294118)+.2e-2*tanh(7.0422535*y[130]-1.3661972)+.155e-1*tanh(77.519380*y[130]-8.1395349)))/ctn
    dy[131] = y[131]-y[85]+y[102]+ (current*Kr)

end
function ocp_cathode(theta_p,up)
    Up = (-5.057-13.6*theta_p^2+121.6*theta_p^4-185.1*theta_p^6-45.43*theta_p^8+127.7*theta_p^10)/(-1-5.04*theta_p^2+32.3*theta_p^4-40.95*theta_p^6-25.94*theta_p^8+40.63*theta_p^10)

end
yyy = function(pars::Array{Float64})
    socn=pars[24]
    socp=pars[25]
    y0 = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10) + 0.5e-1 - 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) + 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) + 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) + 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) + 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) + 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1)]

end
y0=zeros(131)

n=length(y0)
MM=zeros((n,n))
node = 59

for i in 1:node
    MM[i,i]=1.
end

Sorry for all the huge code snippets!

ChrisRackauckas commented 6 years ago

The first run includes compilation which is what the time increase is. I think this is hitting a Julia bug similar to https://github.com/JuliaLang/julia/issues/22255 which is due to an inference recursion limit in v0.6 https://github.com/JuliaLang/julia/issues/26092 . Or it may be due to ForwardDiff.jl compilation time issues like https://github.com/JuliaDiff/ForwardDiff.jl/issues/278 .

To check if it's a ForwardDiff thing, try using Rodas4P(autodiff=false) as the algorithm to instead use numerical differentiation instead of autodifferentiation. If that doesn't help (or maybe even if it does), this might be worth pinging the compiler folks as a real-world example of where compilation time needs some debugging.

nealde commented 6 years ago

Wow, thanks for the super quick response!

Trying it with autodiff=false results in much faster compilation, around 8 seconds, so it definitely seems to be related to that.

ChrisRackauckas commented 6 years ago

I isolated the issue and reported it upstream. Hopefully this gets fixed there.

https://github.com/JuliaLang/julia/issues/27488