MATLAB: Time derivative of parameters within ODE solvers

MATLABode

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.

Best Answer

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
Related Question