A bit of context: I am writing an RK4 ODE solver for Stan, because the existing RK45 and bdf variable-step solvers are incompatible with the ODE I am trying to solve, i.e. :

There is a component (call this comp.A) in the ODE which is dependent on a non-time dependent component (call this comp.B) calculated in the previous time step, as such, comp.B is brought forward by inclusion in the state space (call this X) using an ad-hoc numerical gradient scheme with some step-size h, i.e. pseudo-code:

I see that both RK45(.) and bdf(.) functions takes in functions as arguments, but I cannot find in the documentation how to declare and pass function handles to a different function? Can someone please direct me to the right resource?

This is unfortunately not possible in current version of Stan (RK45 and the like get special treatment from the compiler). It should AFAIK be possible once Stan3 is out, maybe @Bob_Carpenter can clarify?

I was just reading the documentation for Stanâ€™s algebraic solver earlier today and it mentioned that ode solvers and algebraic solver functions receive special treatment with regards to handling functions as inputs, i.e. â€śbasic functionsâ€ť do not have the capability of taking other functions as inputs.

the original intent was to write a general RK4 solver, but I think Iâ€™ll just use a forward declaration and explicitly specify the ODE gradient function in the RK4 solver for now.