Add an option that makes it possible to ensure that the last point is exactly for Tmax; make sure that we never exceed Tmax

ADDENDUM:
`tn` is the beginning of the current integration step `[tn, tn + h]`, so `tn` should start equal to `parameters.initial.time.value`. Moreover, `tn <= parameters.tmax` does not guarantee that we do not overshoot, since it says nothing about `tn + h`.
We need `tn <= parameters.tmax - δt` to guarantee that.