# Issue with BifurcationKit.continuation and arbitrary ODE system - MethodErrors (Convert, ambiguous, Incorrect matrix dimensions...)

**URL:** https://discourse.julialang.org/t/issue-with-bifurcationkit-continuation-and-arbitrary-ode-system-methoderrors-convert-ambiguous-incorrect-matrix-dimensions/107636
**Category:** General Usage
**Tags:** question, error, dynamical-systems
**Created:** [December 14, 2023, 10:44pm UTC](https://discourse.julialang.org/t/issue-with-bifurcationkit-continuation-and-arbitrary-ode-system-methoderrors-convert-ambiguous-incorrect-matrix-dimensions/107636 "2023-12-14T22:44:28Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![Mahnoor\_Pop](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mahnoor_pop/32/205364_2.png) [@Mahnoor\_Pop](https://discourse.julialang.org/u/Mahnoor_Pop)
#### Post date: [December 14, 2023, 10:44pm UTC](https://discourse.julialang.org/t/issue-with-bifurcationkit-continuation-and-arbitrary-ode-system-methoderrors-convert-ambiguous-incorrect-matrix-dimensions/107636/1 "2023-12-14T22:44:28Z")

</div>

Hi all! I’m trying to set up a pipeline for creating an ODE system based on a network of interactions. As a simple example, the corresponding system for the following network would be as such:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/4/f40dab8df5fb20351bb057ec530927ab846bc10d.png)  
\frac{dA}{dt}=-A+\Phi(FB) \\ \frac{dB}{dt}=-B+\Phi(A)  
Where A and B are processes that depend on and amplify each other, and F is some parameter (that I have been trying to bifurcate with). \Phi(x) is a sigmoid function, in my case I use:

```julia
Φ(x, ϵ = 0.05, θ = 0.5) = @. 1.0 / (1.0 + exp(-(x - θ)/ϵ))

```

This is my code for generating my custom ODE equations of motion:

```julia
# create equations of motion
function create_viability_model(odenetworkmodel::ODENetworkModel)
    #should be in place
    return function model!(du, u, p, t=0)
        Φ(x, ϵ = 0.05, θ = 0.5) = @. 1.0 / (1.0 + exp(-(x - θ)/ϵ))
        
        for i in range(1, odenetworkmodel.num_nodes)
            du[i] = -u[i]
            input = 1.0

            for j in range(1, odenetworkmodel.num_nodes)
                if odenetworkmodel.adj_matrix[i,j] == 1
                    input *= u[j]
                end
            end
            
            for k in range(1, length(odenetworkmodel.parameters))
                if odenetworkmodel.par_adj_matrix[k,i] != 0
                    input *= p[k]
                end
            end

            du[i] += Φ(input)
        end

        return du
    end
end

```

ODENetworkModel is a struct I created that holds the information that the user inputs to make such a model. If odenetworkmodel.adj\_matrix[i,j] == 1, then node i depends on node j, so we multiply the input by the current value of node j (I’ve been instructed that if A depends on C, D, E, then the product of the inputs is used for the sigmoid function, i.e. \frac{dA}{dt}=-A+\Phi(CDE). The parameter adjacency matrix has a row for each parameter, so for the kth parameter, if the ith column == 1, then that parameter is being fed into node i.

```julia
struct ODENetworkModel
    num_nodes::Int64
    adj_matrix::Array{Int64}
    par_adj_matrix::Array{Int64}
    parameters::Array{Float64}
    p_index::Int64
    p_start::Float64
    p_end::Float64
    u₀::Array{Float64}
end

```

I have been able to use interactive\_trajectory from DynamicalSystems.jl to visualise the trajectories for when there are 2 nodes (such as the example above), and even find\_attractors(ds) from Attractors.jl seems to be doing fine and finding attractors that make sense for the parameters and context, BUT I can’t seem to get BifurcationKit.continuation to work! I’ve been wrestling with some dimension errors, and am currently stuck on a conversion error:

```julia
function create_bifurcation_diagram(problem)
    p = problem.parameters
    u₀ = problem.u₀

    model! = create_viability_model(problem)

    ds = ContinuousDynamicalSystem(model!, u₀, p)
    # initial_conditions = vec(collect.(Iterators.product(0:0.1:2, 0:0.1:2)))
   
    attractors = find_attractors(ds)
    equilibrium_guess = vec(attractors[1])

!line 76! bfprob = BifurcationProblem(model!, equilibrium_guess, vec(p), (@lens _[problem.p_index]), record_from_solution = recordFromSolution) 

    println(bfprob)
    # continuation options
	opts_br = ContinuationPar(p_min = problem.p_start, p_max = problem.p_end)

    br = BifurcationKit.continuation(bfprob, PALC(), opts_br; bothside = true)

    # scene = Plots.plot(diagram; code = (), title="Viability Bifurcation Diagram", ylim = (0.0,1.0), color=:blue, ylabel = "A")

    # display(scene)
end

function main()
    # don't load from file- does that fix it?
    model = ODENetworkModel(2, [0 1; 1 0], [1 0; 0 0], [0.4, 0.9], 1, 0.1, 0.9, [0.4, 0.5])

    create_bifurcation_diagram(model)
end

Mapping initial conditions to attractors: 100%|███████████████████████████████████████████████████████████████| Time: 0:00:00        
ERROR: MethodError: Cannot `convert` an object of type Float64 to an object of type SVector{2, Float64}

Closest candidates are:
  convert(::Type{SVector{N, T}}, ::CartesianIndex{N}) where {N, T}
   @ StaticArrays C:\Users\smahn\.julia\packages\StaticArrays\yXGNL\src\SVector.jl:8
  convert(::Type{T}, ::Main.LinearAlgebra.Factorization) where T<:AbstractArray
   @ Main.LinearAlgebra c:\Users\smahn\AppData\Local\Programs\Julia-1.9.3\share\julia\stdlib\v1.9\LinearAlgebra\src\factorization.jl:59
  convert(::Type{T}, ::Main.LinearAlgebra.Factorization) where T<:AbstractArray
   @ Main.LinearAlgebra c:\Users\smahn\AppData\Local\Programs\Julia-1.9.3\share\julia\stdlib\v1.9\LinearAlgebra\src\factorization.jl:59
  ...

Stacktrace:
 [1] create_bifurcation_diagram(problem::ODENetworkModel)
   @ Main c:\Users\smahn\JuliaScripts\Week 4\13.12\bifurcation_general.jl:76
 [2] main()
   @ Main c:\Users\smahn\JuliaScripts\Week 4\13.12\bifurcation_general.jl:93
 [3] top-level scope
   @ REPL[5]:1

```

Does anyone have any tips for when you’re defining an arbitrarily large system like this? When I explicitly write out the equations of motion, I can get BifurcationKit.continuation running smoothly, but not with this method.

Many thanks for your help!!