I've searched the internet for weeks for a step in the right direction with this.
I'm trying to solve a first order ODE- not too difficult I know, and I am capable thus far of implementing the Runge Kutta 4th order scheme in general.
However the problem I am having seems to be quite unique due to the nature of the system I am working with.
The ODE I need to integrate is something along the lines of:
dx/dt = f(x(t), y(t), z(t), w(t),...etc)
Where the independent variable (time in this case) does not appear on the RHS of the equation. This is not a problem for me I think. However, each of the time-dependent variables are evaluated by a different code (A stellar astrophysics code) which solves the equations of stellar evolution using an adaptive Runge-Kutta integrator. Hence the output time steps are not all equal. The output from that code is then fed into my Runge-Kutta integrator.
For each iteration in my Runge-Kutta integrator I need to update the time dependent variables at that time, but RK-4 assumes a constant time step. The Runge-Kutta Fehlberg scheme has been suggested to me to solve the problem of varying time steps. But it isn't clear to me how it does in fact solve the problem.
I understand the RKF measures the error at each time step, and compares it to a specified error tolerance and then either increases or decreases the time step accordingly. I have implemented this in my code. My question is: If the RKF is deciding the time step for my integration, yet my input data from the other code still has it's own timestep, how can this really work?!
Every example of either RK-4 or RK Fehlberg I have seen online either assumes input data is evenly spaced or that there is no input data at all. Surely someone has had to solve an ODE where they input data from another solution of ODE's?
The classical RK4, like any other Runge-Kutta method, is a one-step method. Which implies that the step-size has only to be constant within the step and can vary between steps.
The only question then is how to vary the step size. A first variant is to compare one step of size $h$ to two steps of size $h/2$, estimate the local error per Richardson and adapt the step size if that local error deviates too widely from the expected contribution ($tol·h/T$).
Ideally, you would assemble all coupled ODE in one ODE function and solve them simultaneously. If that is not possible or desirable, you provide the values of the other components via interpolation.