Is there a way to implement a pipeline (discrete) delay using ODE solver in desolve package in R?

Question

As part of trying to implement a system dynamics model in R using the DeSolve package, I would like to know if there is a way to build a pipeline (discrete delay) into my model.

An example of a pipeline delay could be a distribution process e.g. where an Amazon package leaves a warehouse at time X and takes 2 days to reach me i.e. the entire package reaches me after a fixed duration of 2 days.

I know that the more popular simulation packages such as Vensim etc. have fixed functions (Delay Fixed etc.) to deal with this but I'm looking for guidance on how this could be implemented in R (including the underlying equation(s) to be able to do so).

If we use a simple example of the kind of thing I'm trying to do - we have a model with:

A. One exogenous variable 1. Desired growth rate = 10%

B. Two stocks:

1. Delay

   • The inflow into this stock is 10% of the value of the Asset stock (see below) in any given period
   • The outflow should be the inflow from two periods ago

2. Asset

   • The inflow into this stock should be the outflow from the Delay stock
   • There is no outflow/decay from this stock

Basically what I'm saying is - if I make an investment in assets at time t=0, this should be realised and reflected in the value of the assets 3 time periods later at the end of t=2.

You'll see in my R code I get stuck at the line where I have to define the delay output equation - mathematically, to achieve the desired delay what I need is to subtract the current delay value from the investment made one time period ago - but I have no idea how to call those lagged values - I've looked into using dede rather than ode but not sure if that's doing what I need it to.

I'm also aware I could achieve a first order delay by dividing the Delay value by 3 in the above equation but this would mean I obtain a third of the benefit of my investment immediately in the current time which is not the effect I'm looking to achieve.

    library(deSolve)

    Start <- 0
    Finish <- 10
    Step <- 1
    simtime <- seq(Start, Finish, by=Step)
    stocks <- c(sAssets=10,sDelay=0)
    auxs <- c(aDesired.Growth=0.10)

    model <- function(time, stocks, auxs) {
    with(as.list(c(stocks,auxs)),
    {

    f.Delay.Input <- sAssets * aDesired.Growth

    f.Delay.Output <- sDelay - [f.Investment #from one time period ago]

    f.Asset.Input <- f.Delay.Output

    da_dt <- f.Asset.Input
    dd_dt <- f.Delay.Input - f.Delay.Output

    return(list(c(da_dt, dd_dt),
    Delay.Input=f.Delay.Input,
    Delay.Output=f.Delay.Output,
    Investment.in.Assets=f.Asset.Input,
    ))
    })
    }

    o <- data.frame(ode(y=stocks, times=simtime, func=model, parms=auxs, 
    method="euler"))

    summary(o)

My expected outcome for the Delay stock outflows are as follows. Note that the inflows are made up since in the actual model, the values will fluctuate in each period:

Time    Delay   Inflow  Outflow 
 `0     0       200     0`
 `1    200      180     0`
 `2    380      80      200 ---> 380 (current delay) -180 (previous inflow)`
 `3    260      176     180 ---> 260-80 `
 `4    256      288     80  ---> 256-176`
 `5    464      XXX     176 ---> 464-288`

..and so on.

Thank you for taking the time to read and if the question has already been asked elsewhere I'd be grateful if you could point me in the right direction as I haven't been able to find an answer.