While solving a system of DAE using the integrator interface, I face a weird behavior with derivatives of the solution vector. At random steps, the resultant value (`da`

here) is returned as non-numeric but a “Dual”. Each time that this happens the solution process freezes. BTW this behavior happens only when autodiff=true.

Anyone got an idea what’s going on?

The first few values of `da`

(dx[2] in the DAE) printed from within `Pᵥ`

whenever it is called are the following:

da=0.0

da=0.0

da=6.9529826154277e-310

da=3.84996053946e-313

da=Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}}(6.649515459556e-312,4.2e-322,6.64930681864e-312,6.64951549331e-312,4.45e-322,6.64930685026e-312,6.649515494497e-312,4.7e-322,6.649306858163e-312,6.64951549568e-312,4.94e-322)

da=0.0

da=0.0

da=0.0

da=0.0

da=0.0

da=0.0

da=Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}}(3.682625513e-315,1.0,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315)

da=-0.06745357092515131

da=-0.27408187162467584

da=-0.41785728598205224

As you can see, sometimes `da`

is non-numeric.

The DAE system is:

```
function TwoPhaseFlow(dx, x, params, t)
p, T, r = params
Tᵢ = T[1]
φ = porosity(x[2], x[3])
rate = -A_g*x[7]*x[8]^(1.22)*exp(-E_g/(R_gas*x[6]))
δ = x[2]/5
dx[1] = (1.0/(ρₛ*x[2]*(1-cbrt(φ))))*(x[8] -p(t) -Pᵥ(x[2], c_radius(α(x[2], x[3]), Tᵢ, p(t)), dx[2], dx[5], x[3], α(x[2], x[3]), T, p(t), r)
-sign(x[1])*Pᵧ(α(x[2], x[3]), T, p(t), r) -2*ρₛ*x[1]*((x[1] + dx[5]/ρₛ)*(1-cbrt(φ))) + 0.5*x[1]^2*(1-cbrt(φ)^4) -dx[5]^2*(1/ρₛ - 1/x[4]))
dx[2] = x[1] + dx[5]/ρₛ
dx[3] = (x[2]/x[3])^2 * (dx[2] - dx[5]/ρₛ)
dx[4] = -3 * (x[4]/x[2]) * (dx[2] - dx[5]/x[4])
dx[5] = (ρₛ * λₛ * Tᵢ * (R_gas*Tᵢ/(E_s)) * A_s * exp(-E_s/(R_gas*Tᵢ))/(q_s*(1-G_unreacted(Tᵢ)) - (Cp_s*(Tᵢ-Tg₀)-q_s)*log(1/G_unreacted(Tᵢ))))^0.5
dx[6] = (1/(x[4]*Cp_mixg)) * ( (3/x[2])*(λg*((Tᵢ - x[6])/δ) + dx[5]*Cp_gR*(Tᵢ-x[6])) - rate*(ΔhfRg - ΔhfP2 + (Cp_gR - Cp_gP2)*(x[6]-Tg₀)) + (1-η_abel*x[4])*dx[8] )
dx[7] = rate/x[4] + (3/(x[2]*x[4]))*dx[5]*(G_unreacted(Tᵢ) - x[7])
dx[8] = x[4]*(R_gas_air)*x[6] / (1-η_abel*x[4]) - x[8]
dx[9] = -rate/x[4] + (3/(x[2]*x[4]))*dx[5]*(1.0 - G_unreacted(Tᵢ) - x[9])
dx[10] = -(3/(x[2]*x[4]))*dx[5]*x[10]
nothing
end
```