Hello, I’d like to ask for your help. In an ODE problem, after I complete the calculations for a given time step ( t ), I need to update a certain state. This state is not time-dependent but is recalculated based on another variable ( u ), and then used for subsequent calculations.

I have a mixed ODE problem described as follows:

```
# where s is in m^3 and q_in, q_out are in m^3/s
ds/dt = q_in - q_out
dq_in = new_q_in - q_in
# where new_q_in is calculated by a function func(q_out), and q_gen is a forcing data.
# We can assume that the function is: func = (q_out, q_gen) -> q_out + q_gen
new_q_in = func(q_out, q_gen)
```

As you can see, these two balance equations are different: one updates the state ( s ) (in m³), and the other updates the flow rate ( q_{\text{in}} ) (in m³/s, which is time-dependent).

I tried to construct an ODE problem based on this formula. The code is as follows:

```
using DataInterpolations
using OrdinaryDiffEq
# q_gen
input_arr = [3, 4, 4, 5, 6, 10, 23, 24, 38, 40, 59, 63, 49, 32, 22, 12, 9, 4]
input_itp = LinearInterpolation(input_arr[1,:], collect(1:length(input_arr)))
# ode function
function discharge!(du, u, p, t)
s_river, q_in = u[1], u[2]
q_out = (s_river + q_in) / (p[1] + 1)
du[1] = q_in - q_out
du[2] = q_out + input_itp(t) - q_in
end
# build problem and solve
prob = DiscreteProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
sol = solve(prob, FunctionMap())
```

In the code, I built a `DiscreteProblem`

to solve this problem. However, mathematically, this formula should be continuous, but when constructing an `ODEProblem`

, the results are different:

```
prob2 = ODEProblem(discharge2!, [0.1, 3], (1, length(input_arr)), [0.3])
sol2 = solve(prob2, Tsit5(), saveat=1.0)
```

The reason lies in the equation `du[2] = q_out + input_itp(t) - q_in`

, which doesn’t hold in the `ODEProblem`

.

Therefore, I attempted to replace ( u[2] ) by using a `DiscreteCallback`

:

```
function condition(u, t, integrator)
return true
end
function affect!(integrator)
q_in = integrator.u[2]
q_out = (integrator.uprev[1] + q_in) / (p[1] + 1)
integrator.u[2] = q_out + input_itp(integrator.t) - q_in
end
cb = DiscreteCallback(condition, affect!)
prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [p])
sol = solve(prob, Tsit5(), dt=1.0, callback=cb)
sol_u = Array(sol)
```

I am not sure if this approach is correct, as the results are still different from earlier (perhaps there’s a mistake in one of the steps).