Hi

I am using julia via R to simulate compartmental epidemic ODE models. Is there a way to return a value that is calculated inside the ODE function for every step of the integration? So far it only returns the states. Below is a toy model, and I would like to return the value of “FOI” and “N” at each time step in addition to the states. What is the best way to do this? `return (du[1], du[2], du[3], FOI, N)`

does not work

```
library(JuliaCall)
library(diffeqr)
func_julia2 <- JuliaCall::julia_eval("
function func_julia2(du, u , p , t)
beta = p[1]
gamma = p[2]
mu = p[3]
sigma = p[4]
eta = p[5]
phi = p[6]
S = u[1]
I = u[2]
R = u[3]
N = S + I + R
beta_eff = beta * (1.0+eta*sin(pi*(t-phi)/365.0)^2.0)
FOI = beta_eff*I/N
du[1] = -FOI*S + mu*N - mu*S + sigma*R
du[2] = FOI*S - gamma*I - mu*I
du[3] = gamma*I - mu*R - sigma*R
end")
JuliaCall::julia_assign("u", c(10000,1,0))
JuliaCall::julia_assign("p", c(0.25, 1/5, 1/(70*365), 1/365, 0.2, 10))
JuliaCall::julia_assign("t", c(0.0, 3650.0))
prob3 = JuliaCall::julia_eval("ODEProblem(func_julia2, u, t, p)")
sol = enviro$solve(prob3, enviro$AutoVern7(enviro$KenCarp4()),
saveat=1.0, abstol=1e-8, reltol=1e-6)
```