# DifferentialEquations.jl - how to calculate an 'assignment' and pass back to an ODE?

**URL:** <https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390>\
**Category:** Modelling & Simulations\
**Tags:** question, diffeq\
**Created:** [October 7, 2021, 7:24pm UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390 "2021-10-07T19:24:41Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![bgsmall](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bgsmall/32/29770_2.png) [@bgsmall](https://discourse.julialang.org/u/bgsmall)\
**Post date:** [October 7, 2021, 7:24pm UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/1 "2021-10-07T19:24:41Z")

</div>

Hi.

I have hit a snag in that I wish to calculate a quantity `u[3]` as an ‘assignment’ (not sure that is the correct terminology) within the ODE function and pass back to one of my ODE’s `du[3]` My code is below:

```julia
# Packages (import Pkg; Pkg.add("name of package")
	using DifferentialEquations
	using Plots 

	function perfadap(du,u,p,t)
	## Reactions
	# Reaction: id = r1
		reaction_r1 = p[1]*p[2]*u[3]; 
	# Reaction: id = r2
		reaction_r2 = p[1]*p[3]*u[1]*u[2];
	# Reaction: id = r3
		reaction_r3 = p[1]*p[4]*u[3]; 
	# Reaction: id = r4
		reaction_r4 = p[1]*p[5]*u[2];

	## Assignment	
	# assignmentRule: variable = S
		u[3] = 1*floor(tspan/p[6]);
	
	## Species	
	# Species: id = R, name = u[1], affected by kineticLaw
		du[1] = (1/(p[1]))*(( 1.0 * reaction_r1) + (-1.0 * reaction_r2));
	# Species: id = X = u, name = X = u[2], affected by kineticLaw
		du[2] = (1/(p[1]))*(( 1.0 * reaction_r3) + (-1.0 * reaction_r4));
	# Species: id = S, name = u[3], involved in a rule 	
	    du[3] = u[3]; 
	end
	
	### PARAMETERS (p)
	p = [1.0, 2.0, 2.0, 1.0, 1.0, 4.0] # compartment vol, k1, k2, k3, k4, tau
	### TIME (tspan)
	tspan = (0.0, 20.0) #10.0
	### INITIAL CONDITIONS (u0)
	u0 = [1.0,0.0,0.0] #du[1] = R, du[2] = X, du[3] = S respectively 
    ### PROB SPECIFICATION (op)	
	op = ODEProblem(perfadap, u0, tspan, p)
	### SYS SOLUTION (sol)
	sol = solve(op, Tsit5()) # use Tsit5 ODE solver
	
	### PLOT
	plot(sol, title = "PerfAdap.jl") # seriestype = :scatter, 

```

This is the full ‘ERROR’ and ‘Stacktrace’ appended below

```julia
ERROR: LoadError: MethodError: no method matching /(::Tuple{Float64, Float64}, ::Float64)
Closest candidates are:
  /(::Any, ::ChainRulesCore.AbstractThunk) at C:\Users\bsmall\.julia\packages\ChainRulesCore\1LqRD\src\differentials\thunks.jl:33
  /(::StridedArray{P, N} where N, ::Real) where P<:Dates.Period at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\Dates\src\deprecated.jl:44
  /(::Union{SparseArrays.SparseVector{Tv, Ti}, SubArray{Tv, 1, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}}, false} where var"#s832"<:SparseArrays.AbstractSparseVector{Tv, Ti}, SubArray{Tv, 1, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false} where var"#s832"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}} where {Tv, Ti}, ::Number) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\SparseArrays\src\sparsevector.jl:1450
  ...
Stacktrace:
  [1] perfadap(du::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64}, t::Float64)
    @ Main z:\Documents\Julia\tutorial\Tyson2003_PerfAdap_V2.jl:23
  [2] ODEFunction
    @ C:\Users\bsmall\.julia\packages\SciMLBase\UIp7W\src\scimlfunctions.jl:334 [inlined]
  [3] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, 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, Vector{Float64}, ODEFunction{true, 
typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), 
Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, 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{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}})
    @ OrdinaryDiffEq C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\perform_step\low_order_rk_perform_step.jl:623
  [4] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\solve.jl:456
  [5] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
    @ OrdinaryDiffEq C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\solve.jl:67
  [6] #__solve#471
    @ C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\solve.jl:4 [inlined]
  [7] __solve
    @ C:\Users\bsmall\.julia\packages\OrdinaryDiffEq\8K0Aj\src\solve.jl:4 [inlined]
  [8] #solve_call#42
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:61 [inlined]
  [9] solve_call
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:48 [inlined]
 [10] #solve_up#44
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:87 [inlined]
 [11] solve_up
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:78 [inlined]
 [12] #solve#43
    @ C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:73 [inlined]
 [13] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, typeof(perfadap), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5)
    @ DiffEqBase C:\Users\bsmall\.julia\packages\DiffEqBase\OPDgm\src\solve.jl:68
 [14] top-level scope
    @ z:\Documents\Julia\tutorial\Tyson2003_PerfAdap_V2.jl:43
in expression starting at z:\Documents\Julia\tutorial\Tyson2003_PerfAdap_V2.jl:43

```

Any and all help / pointers appreciated.

Thanks

Ben

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 7, 2021, 7:48pm UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/2 "2021-10-07T19:48:30Z")

</div>

It is not quite clear to me what you try to do…

```julia
u[3] = 1*floor(tspan/p[6])

```

seems pretty weird – and won’t work. As a minimum, you would have to do `tspan./p[6]` just to make your code work. But it doesn’t seem to make sense… normally, `tspan` should not be available inside `perfadap` – and why would you take the `floor` of a _tuple_ such as `tspan`? Also, you cannot apply the `floor` function to a _tuple_ – you would need to do `floor.(tspan./p[6])` – but why would you do that?

Also, assigning `u[3]` like you do kind of ruins the entire differential equation – at least in a conventional sense.

It would be easier to suggest what to do if you explain what you try to achieve.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 7, 2021, 9:44pm UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/3 "2021-10-07T21:44:44Z")

</div>

> [@BLI](#):
>
> Also, assigning `u[3]` like you do kind of ruins the entire differential equation – at least in a conventional sense.

Yeah, this is not an ODE… I’m confused.

---

<div class="post-metadata">

**Author:** ![bgsmall](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bgsmall/32/29770_2.png) [@bgsmall](https://discourse.julialang.org/u/bgsmall)\
**Post date:** [October 8, 2021, 8:43am UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/4 "2021-10-08T08:43:55Z")

</div>

Dear BLI.

Thanks for your rapid response! The calculation I assign to `u[3]` is a signal (an input) that should trigger the evolution of species `u[1]` and `u[2]` in the context of the model (Both `reaction_r1` and `reaction_r2` are triggered by the output of `u[3]`.

Is this (assignment) what is also called a discontinuity in respect of the system of ODE’s - i.e. the calculation is explicit with respect to time `tspan` but is not explicitly calculated WITHIN the system of ODE’s?

So, in summary what I want to do is calculate the quantity `u[3]` and pass this back to the ODE. I guess if there is a way around it, `u[3]` does not need to remain in the function `perfadap`.

Thanks again for your help, much appreciated!

Ben

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 8, 2021, 11:40am UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/5 "2021-10-08T11:40:27Z")

</div>

But what do you mean by the following expression?

> [@BLI](#):
>
> `u[3] = 1*floor(tspan/p[6])`

What is `tspan` meant to be here? Is it the same `tspan` that you use to define the simulation time span, e.g., `tspan = (0,20.0)`? Or do you mean the _current time_ `t`? Or do you mean the final time `tspan[2] = 20.0`. Or what?

And: why should the value of `u[3]` that triggers the reactions be different from `u[3] = 1*floor(tspan/p[6])`?

---

<div class="post-metadata">

**Author:** ![bgsmall](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bgsmall/32/29770_2.png) [@bgsmall](https://discourse.julialang.org/u/bgsmall)\
**Post date:** [October 8, 2021, 4:29pm UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/6 "2021-10-08T16:29:53Z")

</div>

Hi BLI (Bernt?).

Thanks again for your persistence. `tspan' in this context is the current time (as you correctly point out - `t`).

I have worked on this problem further now and apart from not really understanding how to ‘assign’ the value of a derivative to the expression, i.e. `u[3] = 1*floor(tspan/p[6])` I decided to try a different approach.

To that end, I create a function `sig(t)` that uses nested `ifelse.` to determine when the input is on or off and then re-use this function inside my `perfadap` function.

```julia
# Function that turns signal on and off
        function sig(t) 
            ifelse.(t<4,0,
                ifelse.(t>=4&&t<5,1,
                    ifelse.(t>=8&&t<9,1,
                        ifelse.(t>=12&&t<13,1,
                            ifelse.(t>=16&&t<17,1,0)
                                            )
                                     )
                              )
                       )
        end

```

Output 😉

 ![PerfAdap_V4_R_1](https://global.discourse-cdn.com/julialang/original/3X/9/7/977c4627849bb538c7fec4f9995dcc55369b9c5c.png)

Thanks again (and to Chris R)

Ben

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 8, 2021, 5:45pm UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/7 "2021-10-08T17:45:19Z")

</div>

A couple of additional comments…

1. It is not clear to me how you _used_ your function `sig`, so it is not clear to me how you eventuelly produced the plot.
2. You original model function heading is `function perfadap(du,u,p,t)`. Here, `t` denotes the current time, and you cannot suddenly use `tspan` to denote the current time — so the original code was rather odd.
3. Your reactions are not very clear to me. From your comments, let `u[1]` be the concentration of R, `u[2]` be the concentration of X, and `u[3]` be the concentration (?) of S — what are the “chemical” reactions?

---

<div class="post-metadata">

**Author:** ![vettert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vettert/32/30599_2.png) [@vettert](https://discourse.julialang.org/u/vettert)\
**Post date:** [October 8, 2021, 8:12pm UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/8 "2021-10-08T20:12:24Z")

</div>

Maybe bgsmall is actually looking for the callback functionality of differentialequations.jl?

---

<div class="post-metadata">

**Author:** ![bgsmall](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bgsmall/32/29770_2.png) [@bgsmall](https://discourse.julialang.org/u/bgsmall)\
**Post date:** [October 11, 2021, 9:04am UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/9 "2021-10-11T09:04:28Z")

</div>

Hi BLI.

1. `sig` just replaces instances of `u[3]` in the code.
2. I think (?) that `tspan` represents the time over which I want to see the solution (i.e. the boundaries of the solution to `perfadap`) versus `t` which is current ‘model’ time.
3. The reactions are really just a conceptual model I’m working on.

Thanks again for your help.

Ben

---

<div class="post-metadata">

**Author:** ![bgsmall](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bgsmall/32/29770_2.png) [@bgsmall](https://discourse.julialang.org/u/bgsmall)\
**Post date:** [October 11, 2021, 9:16am UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/10 "2021-10-11T09:16:13Z")

</div>

Hi vettert,

Thank you for your prompt.

When I said ‘assignment’, what I meant in Julia terminology was ‘event’. This example may be of use to others:

[https://diffeq.sciml.ai/stable/features/callback\_functions/#DiscreteCallback-Examples](https://diffeq.sciml.ai/stable/features/callback_functions/#DiscreteCallback-Examples)

I will work on this further.

Thanks

Ben

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 11, 2021, 9:29am UTC](https://discourse.julialang.org/t/differentialequations-jl-how-to-calculate-an-assignment-and-pass-back-to-an-ode/69390/11 "2021-10-11T09:29:22Z")

</div>

> [@bgsmall](#):
>
> `u[3] = 1*floor(tspan/p[6]);`

So – I assume you have removed the above line, and replaced `u[3]` with a call to your `sig` function, i.e., `sig(t)`.

```julia
function perfadap(du,u,p,t)
    # input u computed from sig(t) -- inputs are often denoted u
    u = sig(t)
	## Reactions
	# Reaction: id = r1
		reaction_r1 = p[1]*p[2]*u 
	# Reaction: id = r2
		reaction_r2 = p[1]*p[3]*u[1]*u[2];
	# Reaction: id = r3
		reaction_r3 = p[1]*p[4]*u; 
	# Reaction: id = r4
		reaction_r4 = p[1]*p[5]*u[2];

	## Assignment	
	# assignmentRule: variable = S
	#	u[3] = 1*floor(tspan/p[6]);
	
	## Species	
	# Species: id = R, name = u[1], affected by kineticLaw
		du[1] = (1/(p[1]))*(( 1.0 * reaction_r1) + (-1.0 * reaction_r2));
	# Species: id = X = u, name = X = u[2], affected by kineticLaw
		du[2] = (1/(p[1]))*(( 1.0 * reaction_r3) + (-1.0 * reaction_r4));
	# Species: id = S, name = u[3], involved in a rule 	
	    du[3] = u; 
	end

```

In this case, `sig(t)` is not a state; it’s an input, thus `u = sig(t)`. [Control engineers tend to denote the state by `x`, and not `u`, and inputs by `u` — see above]. Thus, adding a “state” of form `du[3] = u` just integrates up `sig(t)` – which is rarely of interest.

If you had computed the input `u` inside of the function instead of by an external function `sig`, you would probably want to store `u` so that you can view it after the simulation. This can not be done easily in `ODEProblem` formulations, but is straightforward in the slightly more complex form `DAEProblem`.