This is really going to hurt the performance of the ODE solver because it makes the integrand non-smooth. I would parameterize f(t) in a smooth way if possible.
For example, you could parameterize f(t) = \left[\tanh(p(2t/T - 1)) + 1\right] f_\max / 2 where p(x) is a polynomial in a Chebyshev basis on x\in [-1,1] (which allow you to go to very high degree without problems), with unconstrained coefficients. The tanh function ensures that f(t) \in (0, f_\max).
See also ODE non autonomous interpolation needed - #5 by ChrisRackauckas …
Forward-mode AD is also going to hurt you here because its cost scales with the number of parameters, so computing your derivatives here will be ≈ 200x the cost of solving the ODE once. You really want to use reverse-mode (either via AD or via manual adjoint equations — see our course notes, chapter 9 — or with the help of SciMLSensitivity.jl) in this kind of problem if possible.