Imagine some functions like x1[t], y1[t] given as InterpolatingFunction objects, for example as a result of solving an ODE with NDSolve.

Now I have an ODE for x[t], y[t], where (repeated) combinations like x1[t]/Sqrt[(x[t] - x1[t])^2+(y[t] - y1[t])^2] or y1[t]/Sqrt[(x[t] - x1[t])^2+(y[t] - y1[t])^2] appear.

This works fine, however it is inefficient, as for each call of x1[t] or y1[t] the InterpolatingFunction object is evaluated during the NDSolve. What would be much more efficient, is if one could for each time step evaluate x1 and y1 (i.e. each once) and plug in the result. I have not managed so far to make that efficient, as NDSolve does not accept a function of type f[t_?NumericQ, ...] when t is the independent time variable.

Any hints on how to speed up?

Thanks!