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?