# Bifurcation analysis with history in ODEs

**URL:** <https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [July 11, 2022, 7:05pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070 "2022-07-11T19:05:28Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![nkoch](https://avatars.discourse-cdn.com/v4/letter/n/8797f3/32.png) [@nkoch](https://discourse.julialang.org/u/nkoch)\
**Post date:** [July 11, 2022, 7:05pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/1 "2022-07-11T19:05:28Z")

</div>

Hello, I’m relatively new to using Julia and I am trying to do bifurcation analysis.

I have a system of ODEs in which one DE depends on the history of another variable (i.e. du[2] depends on history of u[1]). From what I understand there are a few options how to implement this:

1. Have the history as a parameter and update it using callbacks
2. Have the history as a parameter, use the integrator interface and update it explicitly with each step

I haven’t been able to understand if I can (or how to) use either of these methods with BifurcationKit.jl. Are either of these methods compatible with bifurcation analysis using BifurcationKit.jl? Or is there a better way to implement this?

Thanks!  
Nils

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [July 11, 2022, 10:10pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/2 "2022-07-11T22:10:27Z")

</div>

This sounds like a typical DDE problem

[https://diffeq.sciml.ai/stable/types/dde\_types/](https://diffeq.sciml.ai/stable/types/dde_types/)

Can you say more about the way u\_2 depends on the history of u\_1?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [July 11, 2022, 10:57pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/3 "2022-07-11T22:57:43Z")

</div>

Hi,

It is a DDE problem like said above. For the bifurcation analysis, it depends what you have in mind. If you want to compute the stability of stationary points and detect bifurcations, it would not take me very long. If you want to compute periodic orbits, more work is needed.

---

<div class="post-metadata">

**Author:** ![nkoch](https://avatars.discourse-cdn.com/v4/letter/n/8797f3/32.png) [@nkoch](https://discourse.julialang.org/u/nkoch)\
**Post date:** [July 11, 2022, 11:28pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/4 "2022-07-11T23:28:26Z")

</div>

Thank you both for your responses. The DDE problem looks to be the thing I need. For the history I need a vector of u1 values from t-tau to t. If I understand the documentation correctly I could do this by setting constant\_lags to a vector of lags?

And yes I at the moment I am only interested in the stability of stationary points and the detection of bifurcations.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [July 12, 2022, 5:52am UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/5 "2022-07-12T05:52:23Z")

</div>

Can you post the DDE? some are much easier to handle than others

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [July 12, 2022, 7:28am UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/6 "2022-07-12T07:28:21Z")

</div>

I have a strong suspicion what problem you are trying to solve but it would really help if you could _explicitely_ say how u\_1 depends on u\_2 in your system.

> [@nkoch](#):
>
> For the history I need a vector of u1 values from t-tau to t

This is your attempted solution (which intuitely makes sense). However, instead of telling us your _solution_ it would be more helpful to know more about your _problem_.

See also [https://xyproblem.info/](https://xyproblem.info/)

---

<div class="post-metadata">

**Author:** ![nkoch](https://avatars.discourse-cdn.com/v4/letter/n/8797f3/32.png) [@nkoch](https://discourse.julialang.org/u/nkoch)\
**Post date:** [July 12, 2022, 6:11pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/7 "2022-07-12T18:11:38Z")

</div>

Hi,  
Thank you for your suggestions. I have managed to create a working example of the kind of DDE that I would like to do bifurcation analysis on:

```julia
using Parameters, StatsBase, DifferentialEquations, Plots, DirectionalStatistics, DSP

function FHN_coupled!(du, u, p, t)
    @unpack a, b, c, l, gc = p
    du[1] = u[1] * (a - u[1]) * (u[1] - 1) - u[2] .+ gc.*(u[1] .- u[3])
    du[2] = b * u[1] - c * u[2]
    if t >= l
        du[3] = u[3] * (a - u[3]) * (u[3] - 1) - u[4] .+ gc.*(u[3] .- u[1])
        du[4] = b * u[3] - c * u[4]
        ph1 = hilbert(h[1,:] .- mean(h[1,:]))
        ph2= hilbert(h[2,:] .- mean(h[2,:]))
        u[5] = 1 .- abs.(Circular.var(angle.(ph1) .- angle.(ph2)))
    end
end

u0 = [0.5, 0.0, 0.5, 0.0, 0]
p = (a = -0.1, b = 0.01, c = 0.02, l = 50, gc = 0.025)
tspan = (0.0, 500.0)
datasize = Int(tspan[2])
tsteps = range(tspan[1], tspan[2], length=datasize)
lags = 0:1:200
h = zeros(4, length(lags))
prob = DDEProblem(FHN_coupled!, u0, h, tspan, p; constant_lags=lags)

# run DDE and plot 
sol = solve(prob, Tsit5(), saveat=tsteps)
plot(sol[1,:])
plot!(sol[3,:])
plot!(sol[5,:])

```

 ![FHN_coupled](https://global.discourse-cdn.com/julialang/original/3X/1/3/13b76d40b630899e768e7a4611245d954deeaee8.png)

Would this DDE be easy to handle and to implement?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [July 12, 2022, 9:42pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/8 "2022-07-12T21:42:25Z")

</div>

I dont think this has constant lags

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [July 13, 2022, 12:04pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/9 "2022-07-13T12:04:02Z")

</div>

Please note that your function signature for `FHN_coupled!` does _not_ match the function signature required for a `DDEProblem`.

[https://diffeq.sciml.ai/stable/types/dde\_types/#SciMLBase.DDEProblem](https://diffeq.sciml.ai/stable/types/dde_types/#SciMLBase.DDEProblem)

---

<div class="post-metadata">

**Author:** ![nkoch](https://avatars.discourse-cdn.com/v4/letter/n/8797f3/32.png) [@nkoch](https://discourse.julialang.org/u/nkoch)\
**Post date:** [July 13, 2022, 3:04pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/10 "2022-07-13T15:04:16Z")

</div>

Hi,

What makes you think that this has dependent lags? I don’t think my lags depend on the state of the system.

I have attempted to fix my mistake in my DDEProblem function however I have run into some difficulties:  
I cannot figure out how to pass a vector of lags in to get and array of the history.

```julia
using Parameters, StatsBase, DifferentialEquations, Plots, DirectionalStatistics, DSP

function FHN_coupled(du, u, h, p, t)
    @unpack a, b, c, l, gc, lags = p
    du[1] = u[1] * (a - u[1]) * (u[1] - 1) - u[2] .+ gc.*(u[1] .- u[3])
    du[2] = b * u[1] - c * u[2]
    if t >= l
        du[3] = u[3] * (a - u[3]) * (u[3] - 1) - u[4] .+ gc.*(u[3] .- u[1])
        du[4] = b * u[3] - c * u[4]
        hi = h(p, t .- lags)
        ph1 = hilbert(hi[1,:] .- mean(hi[1,:]))
        ph2= hilbert(hi[3,:] .- mean(hi[3,:]))
        u[5] = 1 .- abs.(Circular.var(angle.(ph1) .- angle.(ph2)))
    end
end

u0 = [0.5, 0.0, 0.5, 0.0, 0]
tspan = (0.0, 500.0)
datasize = Int(tspan[2])
tsteps = range(tspan[1], tspan[2], length=datasize*100)
lags =collect(1.:200.)
p = (a = -0.1, b = 0.01, c = 0.02, l = 50, gc = 0.025, lags=lags)
h(p,t) = zeros(5, length(lags))
tau = [1,2]
prob = DDEProblem(FHN_coupled,u0,h,tspan,p; constant_lags=lags) 

ERROR: LoadError: MethodError: no method matching isless(::Vector{Float64}, ::Float64)
Closest candidates are:
  isless(!Matched::T, ::T) where T<:Union{Float16, Float32, Float64} at C:\Users\nilsk\AppData\Local\Programs\Julia-1.7.2\share\julia\base\float.jl:460
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::AbstractFloat) at C:\Users\nilsk\.julia\packages\StatsBase\pnfAI\src\statmodels.jl:99
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at C:\Users\nilsk\.julia\packages\StatsBase\pnfAI\src\statmodels.jl:90
  ...
Stacktrace:
  [1] <(x::Vector{Float64}, y::Float64)
    @ Base .\operators.jl:352
  [2] (::DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}})(p::NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, t::Vector{Float64}, ::Type{Val{0}}; idxs::Nothing)
    @ DelayDiffEq C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\history_function.jl:26
  [3] (::DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}})(p::NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, t::Vector{Float64}, ::Type{Val{0}}) (repeats 2 times)
    @ DelayDiffEq C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\history_function.jl:21
  [4] FHN_coupled(du::Vector{Float64}, u::Vector{Float64}, h::DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, p::NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, t::Float64)
    @ Main C:\Users\nilsk\AppData\Roaming\JetBrains\PyCharm2021.2\scratches\scratch_52.jl:17
  [5] ODEFunctionWrapper
    @ C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\functionwrapper.jl:59 [inlined]
  [6] perform_step!(integrator::DelayDiffEq.DDEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, DDEProblem{Vector{Float64}, Tuple{Float64, Float64}, Vector{Float64}, Tuple{}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DDEFunction{true, typeof(FHN_coupled), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(h), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DelayDiffEq.DDEStats}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DelayDiffEq.FPSolver{NLFunctional{Rational{Int64}, Rational{Int64}}, true, Float64, DelayDiffEq.FPFunctionalCache{Vector{Float64}, Vector{Float64}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Discontinuity{Float64, Int64}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Float64, Int64, DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Discontinuity{Float64, Int64}}, Vector{Float64}, Float64, Nothing}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, repeat_step::Bool)
    @ OrdinaryDiffEq C:\Users\nilsk\.julia\packages\OrdinaryDiffEq\ZgJ9s\src\perform_step\low_order_rk_perform_step.jl:704
  [7] perform_step!
    @ C:\Users\nilsk\.julia\packages\OrdinaryDiffEq\ZgJ9s\src\perform_step\low_order_rk_perform_step.jl:675 [inlined]
  [8] perform_step!
    @ C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\integrators\interface.jl:125 [inlined]
  [9] solve!(integrator::DelayDiffEq.DDEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, DDEProblem{Vector{Float64}, Tuple{Float64, Float64}, Vector{Float64}, Tuple{}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DDEFunction{true, typeof(FHN_coupled), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(h), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DelayDiffEq.DDEStats}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DelayDiffEq.FPSolver{NLFunctional{Rational{Int64}, Rational{Int64}}, true, Float64, DelayDiffEq.FPFunctionalCache{Vector{Float64}, Vector{Float64}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Discontinuity{Float64, Int64}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Float64, Int64, DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Discontinuity{Float64, Int64}}, Vector{Float64}, Float64, Nothing})
    @ DelayDiffEq C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\solve.jl:358
 [10] #__solve#32
    @ C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\solve.jl:5 [inlined]
 [11] #solve_call#28
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\RHAWf\src\solve.jl:429 [inlined]
 [12] #solve_up#34
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\RHAWf\src\solve.jl:767 [inlined]
 [13] #solve#33
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\RHAWf\src\solve.jl:752 [inlined]
 [14] top-level scope
    @ C:\Users\nilsk\AppData\Roaming\JetBrains\PyCharm2021.2\scratches\scratch_52.jl:36
 [15] include(fname::String)
    @ Base.MainInclude .\client.jl:451
 [16] top-level scope
    @ none:1
in expression starting at C:\Users\nilsk\FHN_coupled_lag.jl:36

```

My current solution to this is to pass the lags as a parameter and to get the history for each lag and construct the array in that manner. I don’t think that this is ideal though, is there a better way to do this?

```julia
using Parameters, StatsBase, DifferentialEquations, Plots, DirectionalStatistics, DSP

function FHN_coupled(du, u, h, p, t)
    @unpack a, b, c, l, gc, lags = p
    du[1] = u[1] * (a - u[1]) * (u[1] - 1) - u[2] .+ gc.*(u[1] .- u[3])
    du[2] = b * u[1] - c * u[2]
    if t >= l
        du[3] = u[3] * (a - u[3]) * (u[3] - 1) - u[4] .+ gc.*(u[3] .- u[1])
        du[4] = b * u[3] - c * u[4]
        hist = zeros(2, length(lags))
        T = typeof(t)
        for i in 1:length(lags)
            hist[1,i] = h(p, t - T(lags[i]))[1]
            hist[2,i] = h(p, t - T(lags[i]))[3]
        end
        ph1 = hilbert(hist[1,:] .- mean(hist[1,:]))
        ph2 = hilbert(hist[2,:] .- mean(hist[2,:]))
        u[5] = 1 .- abs.(Circular.var(angle.(ph1) .- angle.(ph2)))
    end
end

u0 = [0.5, 0.0, 0.5, 0.0, 0]
tspan = (0.0, 500.0)
datasize = Int(tspan[2])
tsteps = range(tspan[1], tspan[2], length=datasize*100)
lags = collect(1:200)
p = (a = -0.1, b = 0.01, c = 0.02, l = 50, gc = 0.025, lags=lags)
h(p,t) = zeros(5)

prob = DDEProblem(FHN_coupled,u0,h,tspan,p; constant_lags=lags)

# run DDE and plot
sol = solve(prob,MethodOfSteps(Tsit5()), saveat=tsteps)
plot(sol[1,:], dpi=300 )
plot!(sol[3,:], dpi=300 )
plot!(sol[5,:], dpi=300 )

```

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [July 13, 2022, 4:04pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/11 "2022-07-13T16:04:51Z")

</div>

My first impression is that

```julia
        for i in 1:length(lags)
            hist[1,i] = h(p, t - T(lags[1]))[1]
            hist[2,i] = h(p, t - T(lags[1]))[3]
        end

```

should probably be

```julia
        for i in 1:length(lags)
            hist[1,i] = h(p, t - T(lags[i]))[1]
            hist[2,i] = h(p, t - T(lags[i]))[3]
        end

```

Is that correct?

Whether this approach is correct in your case, depends on the internals of the function `hilbert`. Could you explain what that does?

Do you really intend to have 200 discrete lags or are you really approximating an integral over the history?

---

<div class="post-metadata">

**Author:** ![nkoch](https://avatars.discourse-cdn.com/v4/letter/n/8797f3/32.png) [@nkoch](https://discourse.julialang.org/u/nkoch)\
**Post date:** [July 13, 2022, 4:32pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/12 "2022-07-13T16:32:52Z")

</div>

Yes, my apologies, that is correct. I have edited the code in my previous reply to reflect this.

As for the function `hilbert(x) ` it returns the analytic signal for x, namely x\_A = x +j \hat{x} where \hat{x} is the Hilbert transform of x. The exact function and its source code can be found here:  
[https://docs.juliadsp.org/stable/util/#DSP.Util.hilbert](https://docs.juliadsp.org/stable/util/#DSP.Util.hilbert)

> <https://github.com/JuliaDSP/DSP.jl/blob/f54471f002551331cf3dbbee3e4e130fddf08204/src/util.jl#L45-L51>

I do intend to have this many lags, as I would like to have a decent estimate of the analytic signal based on the history.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [July 13, 2022, 5:18pm UTC](https://discourse.julialang.org/t/bifurcation-analysis-with-history-in-odes/84070/13 "2022-07-13T17:18:18Z")

</div>

Ok, since you want to use a discrete Hilbert transform, it makes sense to produce discrete samples of the history.

It’s up to you if you want to use quadrature (e.g. QuadGK.jl) to calculate the time average (`mean`) of the history more accurately.