hmax from gsl does not actually enforce a maximum step

about.pl

@@ -8,7 +8,7 @@
     about/0
   ]).
 
-version('4.1.2').
+version('4.1.3').
 
 copyright(
   'Copyright (C) 2003-2018 Inria, EPI Lifeware, Saclay-Île de France, France'

library/gsl_solver.cc

@@ -16,7 +16,7 @@ init_gsl_solver(
     struct gsl_solver *solver, const struct gsl_solver_config *config)
 {
     const gsl_odeiv2_step_type *method = config->method;
-    int variable_count = config->variable_count;
+    size_t variable_count = config->variable_count;
     double *p = (double *) xcalloc(config->parameter_count, sizeof(double));
     gsl_odeiv2_system sys = {functions, jacobian, variable_count, p};
     solver->sys = sys;
@@ -86,7 +86,8 @@ gsl_solver_step(struct gsl_solver_state *state)
     if (t >= config->time_final) {
         return false;
     }
-    t_upper = config->time_final;
+    // hmax not enough to ensure steps smaller than hmax
+    t_upper = std::min(config->time_final, t + solver->d->hmax);
     events(state, &t_upper);
     if (t_upper < t + config->initial_step_size) {
         t_upper = t + config->initial_step_size;

library/gsl_solver.hh

@@ -10,7 +10,7 @@
 
 
 struct gsl_solver_config {
-    int variable_count;
+    size_t variable_count;
     int parameter_count;
     int event_count;
     const gsl_odeiv2_step_type *method;

numerical_simulation.pl

@@ -134,7 +134,7 @@ numerical_simulation :-
   ),
   option(
     maximum_step_size, number, _MaximumStepSize,
-    'maximum step size for the numerical solver, as a fraction of time'
+    'maximum step size for the numerical solver, as a fraction of the total time'
   ),
   option(
     precision, number, _Precision,