DAE/ODE Problem specified in MTK with domains and BCs

[finmod]

@xtalax As you suggested I submit the following DAE problem that I would like to reduce to an ODE system producing a dataset that can be handled by other SciML packages like DataDrivenDiffEq and others.

using ModelingToolkit, DomainSets
using LinearAlgebra
using DataDrivenDiffEq
using OrdinaryDiffEq, Sundials
using Plots; gr()

ModelingToolkit.@parameters begin
    ϕ₀ =0.0401                      # Phillips'curve parameter intercept 0.04/(1. - 0.04^2)
    ϕ₁ =6.41e-05                    # Phillips'curve parameter slope 0.04^3/(1. - 0.04^2)
    η  =0                           # -1 < η < inf 0=CD, inf=Leontief, 1=linear 
    k₀ =-0.0065                     # investment function, constant
    k₁ =-5.0                        # investment function, affine parameter
    k₂ =20                          # investment function exponent
    α  =0.025
    δ  =0.05
    β  =0.02
    r  =0.03
    cor=3.0
    fp =0.333
    b  =0.135
end

@variables begin
    t
    ω(t)  = 0.75
    λ(t)  = 0.90
    d₁(t) = 0.5
    𝛑ₑ(t)
    Y(t)  = 100.0
    𝛎(t)
    𝛟(t) 
    κ(t) 
    g(t)    
end

Dₜ = Differential(t)

domains = [t ∈ Interval(0.0,80.0),
           ω ∈ Interval(0.0,1.0),
           λ ∈ Interval(0.0,1.0),
           d₁ ∈ Interval(0.0,3.0)]


# Variable functions

𝛎 = (t, ω)     -> (1/fp)*((1 -  ω)/b)^η
𝛟 = (t, λ)     -> -ϕ₀ + ϕ₁/((1 - λ)^2)
κ = (t, 𝛑ₑ)    -> k₀ + exp(k₁ + k₂*𝛑ₑ)
g = (t, ω, 𝛑ₑ) -> κ(t, 𝛑ₑ)/𝛎(t, ω) - δ


sfcsys = [    
    Dₜ(ω)  ~ ω*(𝛟(t, λ) - α)
    Dₜ(λ)  ~ λ*(g(t, ω, 𝛑ₑ) - α - β)
    Dₜ(d₁) ~ d₁*(r - κ(t, 𝛑ₑ)/𝛎(t, ω)) + ω - 1 + κ(t, 𝛑ₑ)
    Dₜ(Y)  ~ Y*g(t, ω, 𝛑ₑ)
    0.    ~ 𝛑ₑ - 1 + ω + r*d₁
]

ModelingToolkit.@named primalKeen = ModelingToolkit.ODESystem(sfcsys)
tspan = (0.0, 80.0)

rdKeen = structural_simplify(primalKeen; simplify=true)
full_equations(rdKeen)

dt=0.25
odaeprob = ODAEProblem(rdKeen,[],tspan,jac=true,sparse=true)
odae_sol = solve(odaeprob,Tsit5(),abstol=1/10^14,reltol=1/10^14, saveat=dt)    # 11 ms