Time derivative of parameters within ODE solvers

I have an ODE which has a parameter whose 1st and 2nd order time derivatives are also included in the ODE:

function dy = ODE(t,y)
    f1 = myfun(y(t));
    dy = y + f1 + df1/dt + ddf1/dt^2;
end

Unfortunately, the function of analytical derivatives of myfun is not available. Therefore, df1 and ddf1 can be computed numerically, only. Given that the time step in Matlab ode solvers is not fixed, I wonder if there is a way to numerically compute df1 and ddf1.

Answer 1

[Torsten]

function dy = ODE(t,y)
  dt = 1e-8;
  fm = myfun(t-dt);
  f = myfun(t);
  fp = myfun(t+dt);
  df = (fp - fm) / (2 * dt);
  ddf = (fp - 2 * f + fm) / dt^2;
  dy = y + f + df + ddf;
end

[Seyed Ali Baradaran Birjandi]

Thank you Torsten,

but I am not sure whether we know "dt" using matlab ODE solvers, such as ode45. Because these solvers do not use fixed step time. Please read:

https://www.mathworks.com/matlabcentral/answers/319733-how-to-specifiy-time-in-step-size-using-ode45

[Torsten]

MATLAB's "dt" is irrelevant. As you write, f1 only depends on t and not on y. Thus f1 is independent from the solution variables, and thus its first and second order derivatives at time t can be approximated the way I wrote.

[Seyed Ali Baradaran Birjandi]

Thank you again. My mistake. myfun is a function of y. I will correct it.

Back to the question, is there any way to compute df1 and ddf1?

[Torsten]

Does "myfun" depend explicitly on t or only via y ? I.e. can myfun look like y^2+t^2 or only like y^2 ?

[Seyed Ali Baradaran Birjandi]

It depends on y only, i.e: y^2.

[Torsten]

Let

z = f(y)

the value that "myfun" returns for argument y.

Then

dz/dt = df/dy * dy/dt
d^2z/dt^2 = d^2f/dy^2 * (dy/dt)^2 + df/dy * d^2y/dt^2 

Inserting into your differential equation gives

dy/dt = y + f + df/dy * dy/dt + d^2f/dy^2 * (dy/dt)^2 + df/dy * d^2y/dt^2

or

df/dy * d^2y/dt^2 + (df/dy - 1) * dy/dt + d^2f/dy^2 * (dy/dt)^2 + y + f = 0

Now you can approximate df/dy and d^2f/dy^2 as I suggested above and solve the system (z1 = y, z2 = dy/dt)

z1' = z2
z2' = -((df/dy - 1) * z2 + d^2f/dy^2 * (z2)^2 + z1 + f)/(df/dy)

using ODE45, e.g.