Solving ODE inside a double for loop

General Usage Performance
5

I’m not sure how to properly set up a DAE problem for my case. I have tried the following setup just to see what I get for the kinetic part, but I get an error right at the start.

#display(HTML("<style>.container { width:100% !important; }</style>"))
using Revise
using DifferentialEquations

ρ = 1670.0                                   
Eₐ₁ = 198.74                     
Eₐ₂ = 188.28                      
Eₐ₃ = 143.09                      
lnZ₁ = 45.2                         
lnZ₂ = 40.0                           
lnZ₃ = 34.9                         
R_gas = 0.00831447             


function kinetics(out, dY, Y, p, t)
    TT = p
    R1 = (exp(lnZ₁).*exp(-Eₐ₁/(R_gas.*TT))).*Y[1]       
    R2 = (exp(lnZ₂).*exp(-Eₐ₂/(R_gas.*TT))).*Y[2]       
    R3 = (exp(lnZ₃).*exp(-Eₐ₃/(R_gas.*TT))).*Y[3].^2    

    out[1] = (-1.0/ρ).*R1
    out[2] = (1.0/ρ).*(R1 - R2)
    out[3] = (1.0/ρ).*(R2 - R3)
    out[4] = Y[1] + Y[2] + Y[3] + Y[4] - 1.0
end

u₀ = [1.0, 0.0, 0.0, 0.0]
du₀ = ones(length(u₀))
tspan = (0.0, 1.0e-6)

differential_vars = [true, true, true, false]
chemprob = DAEProblem(kinetics, du₀, u₀, tspan, 1000.0, differential_vars = differential_vars)

using Sundials
sol = solve(chemprob, IDA())

The above code results in the following error:

[IDAS ERROR] IDACalcIC
Newton/Linesearch algorithm failed to converge.
retcode: InitialFailure