Boundary Value Problems in ModelingToolkit

Specific Domains Modelling & Simulations
, ,
1

Is it currently possible to symbolically setup and then numerically solve boundary-value problems using ModelingToolkit? This issue tells me this is possible through the dynamic optimization interface but I’m not sure how to set it up this way. As a test case I’m trying this simple example from the BoundaryValueDiffEq docs:

using BoundaryValueDiffEq
function f!(du, u, p, t)
    du[1] = u[2]
    du[2] = u[1]
end
function bc!(res, u, p, t)
    res[1] = u(0.0)[1] - 1
    res[2] = u(1.0)[1]
end
tspan = (0.0, 1.0)
u0 = [0.0, 0.0]
prob = BVProblem(f!, bc!, u0, tspan)
sol = solve(prob, MIRK4(), dt = 0.01)

Here is my attempt at setting this up in ModelingToolkit:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using BoundaryValueDiffEq # for BVProblem
@variables y(..)

eqs = [D(D(y(t))) ~ y(t)]

(ti, tf) = (0.0, 1.0)
cons = [y(ti) ~ 1.0, y(tf) ~ 0.0]

@named bvpsys = System(eqs, t; constraints=cons)
bvpsys = mtkcompile(bvpsys)

u0 = [0.0, 0.0] 
prob = BVProblem(bvpsys, u0, (ti, tf)) 
# or do I need JuMPDynamicOptProblem() ?
sol = solve(prob, MIRK4(), dt=0.01)

My attempt above leads to a seg fault and my REPL crashes so something is horribly wrong. I’m guessing BVProblem is not the one to use here and I need a dynamic optimization solver like in these examples, but I’m not sure how. What is the proper syntax for this?