I don't think NDSolve has capabilities to use variable boundary limits like that. The only method to try I can think off right now, is to re-scale the time axis to fit this equation. So in your ODEs you introduce tnew = t/tf so that 0 <= tnew <= 1. After that you will have to rethink your ODEs to make sense with that re-scaling and how to get the value for tf in a self-consistent manner.
