diff --git a/include/amici/serialization.h b/include/amici/serialization.h
index b5bfe466fc..97172e8375 100644
--- a/include/amici/serialization.h
+++ b/include/amici/serialization.h
@@ -87,6 +87,7 @@ void serialize(Archive& ar, amici::Solver& s, unsigned int const /*version*/) {
     ar & s.maxtime_;
     ar & s.max_conv_fails_;
     ar & s.max_nonlin_iters_;
+    ar & s.max_step_size_;
 }
 
 /**
diff --git a/include/amici/solver.h b/include/amici/solver.h
index ab4f8ef427..a52e80e82d 100644
--- a/include/amici/solver.h
+++ b/include/amici/solver.h
@@ -964,6 +964,18 @@ class Solver {
      */
     int getMaxConvFails() const;
 
+    /**
+     * @brief Set the maximum step size
+     * @param max_step_size maximum step size. `0.0` means no limit.
+     */
+    void setMaxStepSize(realtype max_step_size);
+
+    /**
+     * @brief Get the maximum step size
+     * @return maximum step size
+     */
+    realtype getMaxStepSize() const;
+
     /**
      * @brief Serialize Solver (see boost::serialization::serialize)
      * @param ar Archive to serialize to
@@ -1752,6 +1764,11 @@ class Solver {
      */
     virtual void apply_max_conv_fails() const = 0;
 
+    /**
+     * @brief Apply the allowed maximum stepsize to the solver.
+     */
+    virtual void apply_max_step_size() const = 0;
+
     /** state (dimension: nx_solver) */
     mutable AmiVector x_{0};
 
@@ -1891,6 +1908,9 @@ class Solver {
      *  step */
     int max_conv_fails_{10};
 
+    /** Maximum allowed step size */
+    realtype max_step_size_{0.0};
+
     /** CPU time, forward solve */
     mutable realtype cpu_time_{0.0};
 
diff --git a/include/amici/solver_cvodes.h b/include/amici/solver_cvodes.h
index 592298b7ec..8abf407f1d 100644
--- a/include/amici/solver_cvodes.h
+++ b/include/amici/solver_cvodes.h
@@ -253,6 +253,8 @@ class CVodeSolver : public Solver {
     void apply_max_nonlin_iters() const override;
 
     void apply_max_conv_fails() const override;
+
+    void apply_max_step_size() const override;
 };
 
 } // namespace amici
diff --git a/include/amici/solver_idas.h b/include/amici/solver_idas.h
index 079cee6399..3d8a9df485 100644
--- a/include/amici/solver_idas.h
+++ b/include/amici/solver_idas.h
@@ -233,6 +233,8 @@ class IDASolver : public Solver {
     void apply_max_nonlin_iters() const override;
 
     void apply_max_conv_fails() const override;
+
+    void apply_max_step_size() const override;
 };
 
 } // namespace amici
diff --git a/src/hdf5.cpp b/src/hdf5.cpp
index 4efed799a8..e980619e5f 100644
--- a/src/hdf5.cpp
+++ b/src/hdf5.cpp
@@ -802,6 +802,11 @@ void writeSolverSettingsToHDF5(
         file.getId(), hdf5Location.c_str(), "maxtime", &dbuffer, 1
     );
 
+    dbuffer = solver.getMaxStepSize();
+    H5LTset_attribute_double(
+        file.getId(), hdf5Location.c_str(), "max_step_size", &dbuffer, 1
+    );
+
     ibuffer = gsl::narrow<int>(solver.getMaxSteps());
     H5LTset_attribute_int(
         file.getId(), hdf5Location.c_str(), "maxsteps", &ibuffer, 1
@@ -1000,6 +1005,12 @@ void readSolverSettingsFromHDF5(
         );
     }
 
+    if (attributeExists(file, datasetPath, "max_step_size")) {
+        solver.setMaxStepSize(
+            getDoubleScalarAttribute(file, datasetPath, "max_step_size")
+        );
+    }
+
     if (attributeExists(file, datasetPath, "maxsteps")) {
         solver.setMaxSteps(getIntScalarAttribute(file, datasetPath, "maxsteps")
         );
diff --git a/src/solver.cpp b/src/solver.cpp
index 88185817dd..352017a6ee 100644
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -46,6 +46,7 @@ Solver::Solver(Solver const& other)
     , check_sensi_steadystate_conv_(other.check_sensi_steadystate_conv_)
     , max_nonlin_iters_(other.max_nonlin_iters_)
     , max_conv_fails_(other.max_conv_fails_)
+    , max_step_size_(other.max_step_size_)
     , maxstepsB_(other.maxstepsB_)
     , sensi_(other.sensi_) {}
 
@@ -195,6 +196,7 @@ void Solver::setup(
 
     apply_max_nonlin_iters();
     apply_max_conv_fails();
+    apply_max_step_size();
 
     cpu_time_ = 0.0;
     cpu_timeB_ = 0.0;
@@ -691,6 +693,14 @@ void Solver::setMaxConvFails(int max_conv_fails) {
 
 int Solver::getMaxConvFails() const { return max_conv_fails_; }
 
+void Solver::setMaxStepSize(realtype max_step_size) {
+    if (max_step_size < 0)
+        throw AmiException("max_step_size must be non-negative.");
+    max_step_size_ = max_step_size;
+}
+
+realtype Solver::getMaxStepSize() const { return max_step_size_; }
+
 int Solver::getNewtonMaxSteps() const { return newton_maxsteps_; }
 
 void Solver::setNewtonMaxSteps(int const newton_maxsteps) {
diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp
index 0952724457..626b4cc050 100644
--- a/src/solver_cvodes.cpp
+++ b/src/solver_cvodes.cpp
@@ -270,6 +270,12 @@ void CVodeSolver::apply_max_conv_fails() const {
         throw CvodeException(status, "CVodeSetMaxConvFails");
 }
 
+void CVodeSolver::apply_max_step_size() const {
+    int status = CVodeSetMaxStep(solver_memory_.get(), getMaxStepSize());
+    if (status != CV_SUCCESS)
+        throw CvodeException(status, "CVodeSetMaxStep");
+}
+
 Solver* CVodeSolver::clone() const { return new CVodeSolver(*this); }
 
 void CVodeSolver::allocateSolver() const {
diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp
index fa08bc29f9..d97d6b3911 100644
--- a/src/solver_idas.cpp
+++ b/src/solver_idas.cpp
@@ -267,6 +267,12 @@ void IDASolver::apply_max_conv_fails() const {
         throw IDAException(status, "IDASetMaxConvFails");
 }
 
+void IDASolver::apply_max_step_size() const {
+    int status = IDASetMaxStep(solver_memory_.get(), getMaxStepSize());
+    if (status != IDA_SUCCESS)
+        throw IDAException(status, "IDASetMaxStep");
+}
+
 Solver* IDASolver::clone() const { return new IDASolver(*this); }
 
 void IDASolver::allocateSolver() const {