I am modelling the human cardiovascular system and solve equations similar to the windkessel equation in the example:

```
# The 5-element WK formulation, also 4-element parallel for Ls = 0.0
function wk5(dP, P, params, t)
#=
The 5-Element WK with serial L
set Ls = 0 for 4-Element WK parallel
Formulation for DifferentialEquations.jl
P: solution vector (pressures p1 and p2)
params: parameter array
[ Rc, Rp, C, Lp, Ls ]
=#
# Split parameter array:
Rc, Rp, C, Lp, Ls = params
dP[1] = (
-Rc / Lp * P[1]
+ (Rc / Lp - 1 / Rp / C) * P[2]
+ Rc * (1 + Ls / Lp) * ForwardDiff.derivative(I, t)
+ I(t) / C
)
dP[2] = -1 / Rp / C * P[2] + I(t) / C
return dP[1], dP[2]
end
```

If I try to solve this, I get an error:

```
┌ Warning: First function call produced NaNs. Exiting.
└ @ OrdinaryDiffEq /home/thor/.julia/packages/OrdinaryDiffEq/T2rXL/src/initdt.jl:81
┌ Warning: Automatic dt set the starting dt as NaN, causing instability.
└ @ OrdinaryDiffEq /home/thor/.julia/packages/OrdinaryDiffEq/T2rXL/src/solve.jl:510
┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
└ @ SciMLBase /home/thor/.julia/packages/SciMLBase/grNUR/src/integrator_interface.jl:325
```

If I use a Central Difference (e.g. the derivative from Calculus.jl), the ODE solver works nicely.

In the above example `I`

is a function that returns a flow rate, for an argument `t`

. `I`

is second order differentiable.

When I compare the analytical, the CD, and the AutoDiff results outside the ODE solver, they are absolutely identical up the 10 decimals.

I have to admit, I am not too familiar with AutoDiff (I do understand the principle mathematically, but not the implementation), so I guess (hope) I’m missing something simple here.