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

`tn` is not overshooting beyond `parameters.tmax`, but the integration step is over `[tn, tn + h]`, so the integration interval overshoots because `tn > parameters.tmax - parameters.Δt`. After that we do `tn += δt`, which is likely greater than `parameters.tmax`, which means we have overshot. Really, the whole confusion comes from the bounds of the loop being wrong.