Imported Upstream version 1.60.0

Change-Id: Ie709530d6d5841088ceaba025cbe175a4ef43050
Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>

28 files changed, 1074 insertions, 271 deletions

diff --git a/boost/numeric/odeint/algebra/default_operations.hpp b/boost/numeric/odeint/algebra/default_operations.hpp
index b10944fa19..56139083cc 100644
--- a/boost/numeric/odeint/algebra/default_operations.hpp
+++ b/boost/numeric/odeint/algebra/default_operations.hpp
@@ -529,7 +529,7 @@ struct default_operations
         { }
 
         template< class Res , class T1 , class T2 , class T3 , class T4 >
-        Res operator()( Res r , const T1 &x_old , const T2 &x , const T3 &dxdt_old , const T4 &x_err )
+        Res operator()( Res r , const T1 &x_old , const T2 &/*x*/ , const T3 &dxdt_old , const T4 &x_err )
         {
             BOOST_USING_STD_MAX();
             using std::abs;
@@ -574,7 +574,7 @@ struct default_operations
         { }
 
         template< class Res , class T1 , class T2 , class T3 , class T4 >
-        Res operator()( Res r , const T1 &x_old , const T2 &x , const T3 &dxdt_old , const T4 &x_err )
+        Res operator()( Res r , const T1 &x_old , const T2 &/*x*/ , const T3 &dxdt_old , const T4 &x_err )
         {
             using std::abs;
             Res tmp = abs( get_unit_value( x_err ) ) /
diff --git a/boost/numeric/odeint/integrate/check_adapter.hpp b/boost/numeric/odeint/integrate/check_adapter.hpp
new file mode 100644
index 0000000000..3d3ebd6c88
--- /dev/null
+++ b/boost/numeric/odeint/integrate/check_adapter.hpp
@@ -0,0 +1,222 @@
+/*
+ [auto_generated]
+ boost/numeric/odeint/integrate/check_adapter.hpp
+
+ [begin_description]
+ Adapters to add checking facility to stepper and observer
+ [end_description]
+
+ Copyright 2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+#ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_CHECK_ADAPTER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_INTEGRATE_CHECK_ADAPTER_HPP_INCLUDED
+
+#include <boost/numeric/odeint/stepper/stepper_categories.hpp>
+#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+template<class Stepper, class Checker,
+         class StepperCategory = typename base_tag<typename Stepper::stepper_category>::type>
+class checked_stepper;
+
+
+/**
+ * \brief Adapter to combine basic stepper and checker.
+ */
+template<class Stepper, class Checker>
+class checked_stepper<Stepper, Checker, stepper_tag>
+{
+
+public:
+    typedef Stepper stepper_type;
+    typedef Checker checker_type;
+    // forward stepper typedefs
+    typedef typename stepper_type::state_type state_type;
+    typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::deriv_type deriv_type;
+    typedef typename stepper_type::time_type time_type;
+
+private:
+    stepper_type &m_stepper;
+    checker_type &m_checker;
+
+public:
+    /**
+     * \brief Construct the checked_stepper.
+     */
+    checked_stepper(stepper_type &stepper, checker_type &checker)
+            : m_stepper(stepper), m_checker(checker) { }
+
+    /**
+     * \brief forward of the do_step method
+     */
+    template<class System, class StateInOut>
+    void do_step(System system, StateInOut &state, const time_type t, const time_type dt)
+    {
+        // do the step
+        m_stepper.do_step(system, state, t, dt);
+        // call the checker
+        m_checker();
+    }
+};
+
+
+/**
+ * \brief Adapter to combine controlled stepper and checker.
+ */
+template<class ControlledStepper, class Checker>
+class checked_stepper<ControlledStepper, Checker, controlled_stepper_tag>
+{
+
+public:
+    typedef ControlledStepper stepper_type;
+    typedef Checker checker_type;
+    // forward stepper typedefs
+    typedef typename stepper_type::state_type state_type;
+    typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::deriv_type deriv_type;
+    typedef typename stepper_type::time_type time_type;
+
+private:
+    stepper_type &m_stepper;
+    checker_type &m_checker;
+
+public:
+    /**
+     * \brief Construct the checked_stepper.
+     */
+    checked_stepper(stepper_type &stepper, checker_type &checker)
+            : m_stepper(stepper), m_checker(checker) { }
+
+    /**
+     * \brief forward of the do_step method
+     */
+    template< class System , class StateInOut >
+    controlled_step_result try_step( System system , StateInOut &state , time_type &t , time_type &dt )
+    {
+        // do the step
+        if( m_stepper.try_step(system, state, t, dt) == success )
+        {
+            // call the checker if step was successful
+            m_checker();
+            return success;
+        } else
+        {
+            // step failed -> return fail
+            return fail;
+        }
+    }
+};
+
+
+/**
+ * \brief Adapter to combine dense out stepper and checker.
+ */
+template<class DenseOutStepper, class Checker>
+class checked_stepper<DenseOutStepper, Checker, dense_output_stepper_tag>
+{
+
+public:
+    typedef DenseOutStepper stepper_type;
+    typedef Checker checker_type;
+    // forward stepper typedefs
+    typedef typename stepper_type::state_type state_type;
+    typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::deriv_type deriv_type;
+    typedef typename stepper_type::time_type time_type;
+
+private:
+    stepper_type &m_stepper;
+    checker_type &m_checker;
+
+public:
+    /**
+     * \brief Construct the checked_stepper.
+     */
+    checked_stepper(stepper_type &stepper, checker_type &checker)
+            : m_stepper(stepper), m_checker(checker) { }
+
+
+    template< class System >
+    std::pair< time_type , time_type > do_step( System system )
+    {
+        m_checker();
+        return m_stepper.do_step(system);
+    }
+
+    /* provide the remaining dense out stepper interface */
+    template< class StateType >
+    void initialize( const StateType &x0 , time_type t0 , time_type dt0 )
+    { m_stepper.initialize(x0, t0, dt0); }
+
+
+    template< class StateOut >
+    void calc_state( time_type t , StateOut &x ) const
+    { m_stepper.calc_state(t, x); }
+
+    template< class StateOut >
+    void calc_state( time_type t , const StateOut &x ) const
+    { m_stepper.calc_state(t, x); }
+
+    const state_type& current_state( void ) const
+    { return m_stepper.current_state(); }
+
+    time_type current_time( void ) const
+    { return m_stepper.current_time(); }
+
+    const state_type& previous_state( void ) const
+    { return m_stepper.previous_state(); }
+
+    time_type previous_time( void ) const
+    { return m_stepper.previous_time(); }
+
+    time_type current_time_step( void ) const
+    { return m_stepper.current_time_step(); }
+
+};
+
+
+/**
+ * \brief Adapter to combine observer and checker.
+ */
+template<class Observer, class Checker>
+class checked_observer
+{
+public:
+    typedef Observer observer_type;
+    typedef Checker checker_type;
+
+private:
+    observer_type &m_observer;
+    checker_type &m_checker;
+
+public:
+    checked_observer(observer_type &observer, checker_type &checker)
+            : m_observer(observer), m_checker(checker)
+    {}
+
+    template< class State , class Time >
+    void operator()(const State& state, Time t) const
+    {
+        // call the observer
+        m_observer(state, t);
+        // reset the checker
+        m_checker.reset();
+    }
+};
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+#endif
\ No newline at end of file
diff --git a/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp b/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp
index 743e57709c..7516d44087 100644
--- a/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp
+++ b/boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp
@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
 
  Distributed under the Boost Software License, Version 1.0.
@@ -25,6 +25,7 @@
 
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
 #include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
 #include <boost/numeric/odeint/util/bind.hpp>
 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
@@ -41,7 +42,7 @@ namespace odeint {
 namespace detail {
 
 // forward declaration
-template< class Stepper , class System , class State , class Time , class Observer>
+template< class Stepper , class System , class State , class Time , class Observer >
 size_t integrate_const(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time end_time , Time dt ,
@@ -74,7 +75,7 @@ size_t integrate_adaptive(
 
 
 /*
- * classical integrate adaptive
+ * integrate adaptive for controlled stepper
  */
 template< class Stepper , class System , class State , class Time , class Observer >
 size_t integrate_adaptive(
@@ -86,8 +87,7 @@ size_t integrate_adaptive(
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
 
-    const size_t max_attempts = 1000;
-    const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
+    failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
     size_t count = 0;
     while( less_with_sign( start_time , end_time , dt ) )
     {
@@ -97,15 +97,14 @@ size_t integrate_adaptive(
             dt = end_time - start_time;
         }
 
-        size_t trials = 0;
         controlled_step_result res;
         do
         {
             res = st.try_step( system , start_state , start_time , dt );
-            ++trials;
+            fail_checker();  // check number of failed steps
         }
-        while( ( res == fail ) && ( trials < max_attempts ) );
-        if( trials == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ) );
+        while( res == fail );
+        fail_checker.reset();  // if we reach here, the step was successful -> reset fail checker
 
         ++count;
     }
diff --git a/boost/numeric/odeint/integrate/detail/integrate_const.hpp b/boost/numeric/odeint/integrate/detail/integrate_const.hpp
index 312564ff16..7a86b32fa6 100644
--- a/boost/numeric/odeint/integrate/detail/integrate_const.hpp
+++ b/boost/numeric/odeint/integrate/detail/integrate_const.hpp
@@ -6,7 +6,7 @@
  integrate const implementation
  [end_description]
 
- Copyright 2012 Mario Mulansky
+ Copyright 2012-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
  Copyright 2012 Karsten Ahnert
 
@@ -43,12 +43,13 @@ template< class Stepper , class System , class State , class Time , class Observ
 size_t integrate_const(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time end_time , Time dt ,
-        Observer observer , stepper_tag 
+        Observer observer , stepper_tag
 )
 {
+
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
-    
+
     Time time = start_time;
     int step = 0;
     // cast time+dt explicitely in case of expression templates (e.g. multiprecision)
@@ -72,11 +73,11 @@ template< class Stepper , class System , class State , class Time , class Observ
 size_t integrate_const(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time end_time , Time dt ,
-        Observer observer , controlled_stepper_tag 
+        Observer observer , controlled_stepper_tag
 )
 {
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
-    
+
     Time time = start_time;
     const Time time_step = dt;
     int real_steps = 0;
@@ -85,8 +86,10 @@ size_t integrate_const(
     while( less_eq_with_sign( static_cast<Time>(time+time_step) , end_time , dt ) )
     {
         obs( start_state , time );
-        real_steps += detail::integrate_adaptive( stepper , system , start_state , time , time+time_step , dt ,
-                                    null_observer() , controlled_stepper_tag() );
+        // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
+        real_steps += detail::integrate_adaptive(stepper, system, start_state, time,
+                                                 static_cast<Time>(time + time_step), dt,
+                                                 null_observer(), controlled_stepper_tag());
         // direct computation of the time avoids error propagation happening when using time += dt
         // we need clumsy type analysis to get boost units working here
         step++;
@@ -102,12 +105,12 @@ template< class Stepper , class System , class State , class Time , class Observ
 size_t integrate_const(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time end_time , Time dt ,
-        Observer observer , dense_output_stepper_tag 
+        Observer observer , dense_output_stepper_tag
 )
 {
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
-    
+
     Time time = start_time;
     
     st.initialize( start_state , time , dt );
@@ -117,7 +120,7 @@ size_t integrate_const(
     int obs_step( 1 );
     int real_step( 0 );
     
-    while( less_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
+    while( less_eq_with_sign( static_cast<Time>(time+dt) , end_time , dt ) )
     {
         while( less_eq_with_sign( time , st.current_time() , dt ) )
         {
@@ -148,6 +151,7 @@ size_t integrate_const(
         
     }
     // last observation, if we are still in observation interval
+    // might happen due to finite precision problems when computing the the time
     if( less_eq_with_sign( time , end_time , dt ) )
     {
         st.calc_state( time , start_state );
diff --git a/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp b/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp
index 3c1d171620..2ef490d592 100644
--- a/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp
+++ b/boost/numeric/odeint/integrate/detail/integrate_n_steps.hpp
@@ -6,7 +6,7 @@
  integrate steps implementation
  [end_description]
 
- Copyright 2012 Mario Mulansky
+ Copyright 2012-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
  Copyright 2012 Karsten Ahnert
 
@@ -32,10 +32,10 @@ namespace detail {
 
 // forward declaration
 template< class Stepper , class System , class State , class Time , class Observer >
-size_t integrate_adaptive(
+size_t integrate_adaptive_checked(
         Stepper stepper , System system , State &start_state ,
         Time &start_time , Time end_time , Time &dt ,
-        Observer observer , controlled_stepper_tag
+        Observer observer, controlled_stepper_tag
 );
 
 
@@ -66,7 +66,7 @@ Time integrate_n_steps(
 
 
 /* controlled version */
-template< class Stepper , class System , class State , class Time , class Observer>
+template< class Stepper , class System , class State , class Time , class Observer >
 Time integrate_n_steps(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time dt , size_t num_of_steps ,
@@ -80,8 +80,9 @@ Time integrate_n_steps(
     for( size_t step = 0; step < num_of_steps ; ++step )
     {
         obs( start_state , time );
-        detail::integrate_adaptive( stepper , system , start_state , time , static_cast<Time>(time+time_step) , dt ,
-                null_observer() , controlled_stepper_tag() );
+        // integrate_adaptive_checked uses the given checker to throw if an overflow occurs
+        detail::integrate_adaptive(stepper, system, start_state, time, static_cast<Time>(time + time_step), dt,
+                                   null_observer(), controlled_stepper_tag());
         // direct computation of the time avoids error propagation happening when using time += dt
         // we need clumsy type analysis to get boost units working here
         time = start_time + static_cast< typename unit_value_type<Time>::type >(step+1) * time_step;
@@ -93,7 +94,7 @@ Time integrate_n_steps(
 
 
 /* dense output version */
-template< class Stepper , class System , class State , class Time , class Observer>
+template< class Stepper , class System , class State , class Time , class Observer >
 Time integrate_n_steps(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time dt , size_t num_of_steps ,
@@ -135,6 +136,7 @@ Time integrate_n_steps(
         }
     }
 
+    // make sure we really end exactly where we should end
     while( st.current_time() < end_time )
     {
         if( less_with_sign( end_time ,
@@ -144,7 +146,7 @@ Time integrate_n_steps(
         st.do_step( system );
     }
 
-    // observation at end point, only if we ended exactly on the end-point (or above due to finite precision)
+    // observation at final point
     obs( st.current_state() , end_time );
 
     return time;
diff --git a/boost/numeric/odeint/integrate/detail/integrate_times.hpp b/boost/numeric/odeint/integrate/detail/integrate_times.hpp
index d5446ba590..2e27990412 100644
--- a/boost/numeric/odeint/integrate/detail/integrate_times.hpp
+++ b/boost/numeric/odeint/integrate/detail/integrate_times.hpp
@@ -6,7 +6,7 @@
  Default integrate times implementation.
  [end_description]
 
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2012 Karsten Ahnert
  Copyright 2012 Christoph Koke
 
@@ -26,6 +26,7 @@
 #include <boost/numeric/odeint/util/unwrap_reference.hpp>
 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
 
 
 namespace boost {
@@ -38,15 +39,18 @@ namespace detail {
 /*
  * integrate_times for simple stepper
  */
-template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
+template<class Stepper, class System, class State, class TimeIterator, class Time, class Observer>
 size_t integrate_times(
         Stepper stepper , System system , State &start_state ,
         TimeIterator start_time , TimeIterator end_time , Time dt ,
         Observer observer , stepper_tag
 )
 {
-    typename odeint::unwrap_reference< Observer >::type &obs = observer;
-    typename odeint::unwrap_reference< Stepper >::type &st = stepper;
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+
+    stepper_type &st = stepper;
+    observer_type &obs = observer;
     typedef typename unit_value_type<Time>::type time_type;
 
     size_t steps = 0;
@@ -82,12 +86,10 @@ size_t integrate_times(
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
     typedef typename unit_value_type<Time>::type time_type;
 
-    const size_t max_attempts = 1000;
-    const char *error_string = "Integrate adaptive : Maximal number of iterations reached. A step size could not be found.";
+    failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
     size_t steps = 0;
     while( true )
     {
-        size_t fail_steps = 0;
         Time current_time = *start_time++;
         obs( start_state , current_time );
         if( start_time == end_time )
@@ -99,15 +101,16 @@ size_t integrate_times(
             if( st.try_step( system , start_state , current_time , current_dt ) == success )
             {
                 ++steps;
+                // successful step -> reset the fail counter, see #173
+                fail_checker.reset();
                 // continue with the original step size if dt was reduced due to observation
                 dt = max_abs( dt , current_dt );
             }
             else
             {
-                ++fail_steps;
+                fail_checker();  // check for possible overflow of failed steps in step size adjustment
                 dt = current_dt;
             }
-            if( fail_steps == max_attempts ) BOOST_THROW_EXCEPTION( std::overflow_error( error_string ));
         }
     }
     return steps;
@@ -125,6 +128,7 @@ size_t integrate_times(
 {
     typename odeint::unwrap_reference< Observer >::type &obs = observer;
     typename odeint::unwrap_reference< Stepper >::type &st = stepper;
+
     typedef typename unit_value_type<Time>::type time_type;
 
     if( start_time == end_time )
diff --git a/boost/numeric/odeint/integrate/integrate_adaptive.hpp b/boost/numeric/odeint/integrate/integrate_adaptive.hpp
index 12bdc8237d..09997142fc 100644
--- a/boost/numeric/odeint/integrate/integrate_adaptive.hpp
+++ b/boost/numeric/odeint/integrate/integrate_adaptive.hpp
@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
 
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or
diff --git a/boost/numeric/odeint/integrate/integrate_const.hpp b/boost/numeric/odeint/integrate/integrate_const.hpp
index ba2d891e97..fae683f978 100644
--- a/boost/numeric/odeint/integrate/integrate_const.hpp
+++ b/boost/numeric/odeint/integrate/integrate_const.hpp
@@ -8,7 +8,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
 
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or
@@ -23,6 +23,7 @@
 
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 #include <boost/numeric/odeint/integrate/null_observer.hpp>
+#include <boost/numeric/odeint/integrate/check_adapter.hpp>
 #include <boost/numeric/odeint/integrate/detail/integrate_const.hpp>
 #include <boost/numeric/odeint/integrate/detail/integrate_adaptive.hpp>
 
@@ -31,91 +32,120 @@ namespace numeric {
 namespace odeint {
 
 
-
-
-
 /*
  * Integrates with constant time step dt.
  */
-template< class Stepper , class System , class State , class Time , class Observer >
+template<class Stepper, class System, class State, class Time, class Observer, class StepOverflowChecker>
 size_t integrate_const(
-        Stepper stepper , System system , State &start_state ,
-        Time start_time , Time end_time , Time dt ,
-        Observer observer
-)
-{
-    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+        Stepper stepper, System system, State &start_state,
+        Time start_time, Time end_time, Time dt,
+        Observer observer, StepOverflowChecker checker
+) {
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
     // we want to get as fast as possible to the end
-    if( boost::is_same< null_observer , Observer >::value )
-    {
+    // no overflow checks needed
+    if (boost::is_same<null_observer, Observer>::value) {
         return detail::integrate_adaptive(
-                stepper , system , start_state ,
-                start_time , end_time  , dt ,
-                observer , stepper_category() );
+                stepper, system, start_state,
+                start_time, end_time, dt,
+                observer, stepper_category());
+    }
+    else {
+        // unwrap references
+        typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+        typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+        typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+
+        return detail::integrate_const(checked_stepper<stepper_type, checker_type>(stepper, checker),
+                                       system, start_state,
+                                       start_time, end_time, dt,
+                                       checked_observer<observer_type, checker_type>(observer, checker),
+                                       stepper_category());
     }
-    else
-    {
-        return detail::integrate_const( stepper , system , start_state , 
-                                        start_time , end_time , dt ,
-                                        observer , stepper_category() );
-      }
 }
 
 /**
- * \brief Second version to solve the forwarding problem, 
- * can be called with Boost.Range as start_state.
- */
-template< class Stepper , class System , class State , class Time , class Observer >
+* \brief Second version to solve the forwarding problem,
+* can be called with Boost.Range as start_state.
+*/
+template<class Stepper, class System, class State, class Time, class Observer, class StepOverflowChecker >
 size_t integrate_const(
-        Stepper stepper , System system , const State &start_state ,
-        Time start_time , Time end_time , Time dt ,
-        Observer observer
-)
-{
-    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+        Stepper stepper, System system, const State &start_state,
+        Time start_time, Time end_time, Time dt,
+        Observer observer, StepOverflowChecker checker
+) {
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
     // we want to get as fast as possible to the end
-    if( boost::is_same< null_observer , Observer >::value )
-    {
+
+    if (boost::is_same<null_observer, Observer>::value) {
         return detail::integrate_adaptive(
-                stepper , system , start_state ,
-                start_time , end_time  , dt ,
-                observer , stepper_category() );
+                stepper, system, start_state,
+                start_time, end_time, dt,
+                observer, stepper_category());
     }
-    else
-    {
-        return detail::integrate_const( stepper , system , start_state , 
-                                        start_time , end_time , dt ,
-                                        observer , stepper_category() );
+    else {
+        typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+        typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+        typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+
+        return detail::integrate_const(checked_stepper<stepper_type, checker_type>(stepper, checker),
+                                       system, start_state,
+                                       start_time, end_time, dt,
+                                       checked_observer<observer_type, checker_type>(observer, checker),
+                                       stepper_category());
     }
 }
 
 
+/**
+* \brief integrate_const without step overflow checker
+*/
+template<class Stepper, class System, class State, class Time, class Observer>
+size_t integrate_const(
+        Stepper stepper, System system, State &start_state,
+        Time start_time, Time end_time, Time dt, Observer observer)
+{
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
+    return detail::integrate_const(stepper, system, start_state,
+                                   start_time, end_time, dt, observer, stepper_category());
+}
 
+/**
+* \brief Second version to solve the forwarding problem,
+* can be called with Boost.Range as start_state.
+*/
+template<class Stepper, class System, class State, class Time, class Observer>
+size_t integrate_const(
+        Stepper stepper, System system, const State &start_state,
+        Time start_time, Time end_time, Time dt, Observer observer
+) {
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
+    return detail::integrate_const(stepper, system, start_state,
+                                   start_time, end_time, dt, observer, stepper_category());
+}
 
 
 /**
- * \brief integrate_const without observer calls
- */
-template< class Stepper , class System , class State , class Time >
+* \brief integrate_const without observer calls
+*/
+template<class Stepper, class System, class State, class Time>
 size_t integrate_const(
-        Stepper stepper , System system , State &start_state ,
-        Time start_time , Time end_time , Time dt
-)
-{
-    return integrate_const( stepper , system , start_state , start_time , end_time , dt , null_observer() );
+        Stepper stepper, System system, State &start_state,
+        Time start_time, Time end_time, Time dt
+) {
+    return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer());
 }
 
 /**
- * \brief Second version to solve the forwarding problem,
- * can be called with Boost.Range as start_state.
- */
-template< class Stepper , class System , class State , class Time >
+* \brief Second version to solve the forwarding problem,
+* can be called with Boost.Range as start_state.
+*/
+template<class Stepper, class System, class State, class Time>
 size_t integrate_const(
-        Stepper stepper , System system , const State &start_state ,
-        Time start_time , Time end_time , Time dt
-)
-{
-    return integrate_const( stepper , system , start_state , start_time , end_time , dt , null_observer() );
+        Stepper stepper, System system, const State &start_state,
+        Time start_time, Time end_time, Time dt
+) {
+    return integrate_const(stepper, system, start_state, start_time, end_time, dt, null_observer());
 }
 
 
@@ -124,30 +154,37 @@ size_t integrate_const(
 
 
 /********* DOXYGEN *********/
-    /**
-     * \fn integrate_const( Stepper stepper , System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
-     * \brief Integrates the ODE with constant step size.
-     *
-     * Integrates the ODE defined by system using the given stepper.
-     * This method ensures that the observer is called at constant intervals dt.
-     * If the Stepper is a normal stepper without step size control, dt is also
-     * used for the numerical scheme. If a ControlledStepper is provided, the 
-     * algorithm might reduce the step size to meet the error bounds, but it is 
-     * ensured that the observer is always called at equidistant time points
-     * t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
-     * and the dense output is used to call the observer at equidistant time
-     * points.
-     *
-     * \param stepper The stepper to be used for numerical integration.
-     * \param system Function/Functor defining the rhs of the ODE.
-     * \param start_state The initial condition x0.
-     * \param start_time The initial time t0.
-     * \param end_time The final integration time tend.
-     * \param dt The time step between observer calls, _not_ necessarily the 
-     * time step of the integration.
-     * \param observer Function/Functor called at equidistant time intervals.
-     * \return The number of steps performed.
-     */
+/**
+ * \fn integrate_const( Stepper stepper , System system , State &start_state , Time start_time ,
+ *                      Time end_time , Time dt , Observer observer , StepOverflowChecker checker )
+ * \brief Integrates the ODE with constant step size.
+ *
+ * Integrates the ODE defined by system using the given stepper.
+ * This method ensures that the observer is called at constant intervals dt.
+ * If the Stepper is a normal stepper without step size control, dt is also
+ * used for the numerical scheme. If a ControlledStepper is provided, the
+ * algorithm might reduce the step size to meet the error bounds, but it is
+ * ensured that the observer is always called at equidistant time points
+ * t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
+ * and the dense output is used to call the observer at equidistant time
+ * points.
+ * If a max_step_checker is provided as StepOverflowChecker, a
+ * no_progress_error is thrown if too many steps (default: 500) are performed
+ * without progress, i.e. in between observer calls. If no checker is provided,
+ * no such overflow check is performed.
+ *
+ * \param stepper The stepper to be used for numerical integration.
+ * \param system Function/Functor defining the rhs of the ODE.
+ * \param start_state The initial condition x0.
+ * \param start_time The initial time t0.
+ * \param end_time The final integration time tend.
+ * \param dt The time step between observer calls, _not_ necessarily the
+ * time step of the integration.
+ * \param observer [optional] Function/Functor called at equidistant time intervals.
+ * \param checker [optional] Functor to check for step count overflows, if no
+ * checker is provided, no exception is thrown.
+ * \return The number of steps performed.
+ */
 
 } // namespace odeint
 } // namespace numeric
diff --git a/boost/numeric/odeint/integrate/integrate_n_steps.hpp b/boost/numeric/odeint/integrate/integrate_n_steps.hpp
index f85c453bd2..7f3a49bddc 100644
--- a/boost/numeric/odeint/integrate/integrate_n_steps.hpp
+++ b/boost/numeric/odeint/integrate/integrate_n_steps.hpp
@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
 
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or
@@ -34,35 +34,80 @@ namespace odeint {
  *
  * the two overloads are needed in order to solve the forwarding problem
  */
-template< class Stepper , class System , class State , class Time , class Observer>
+template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
 Time integrate_n_steps(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time dt , size_t num_of_steps ,
-        Observer observer )
+        Observer observer , StepOverflowChecker checker )
 {
-    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+    // unwrap references
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+    typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+    typedef typename stepper_type::stepper_category stepper_category;
+
     return detail::integrate_n_steps(
-                stepper , system , start_state ,
+                checked_stepper<stepper_type, checker_type>(stepper, checker),
+                system , start_state ,
                 start_time , dt , num_of_steps ,
-                observer , stepper_category() );
+                checked_observer<observer_type, checker_type>(observer, checker),
+                stepper_category() );
 }
 
 /**
  * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
  */
-template< class Stepper , class System , class State , class Time , class Observer >
+template< class Stepper , class System , class State , class Time , class Observer , class StepOverflowChecker >
 Time integrate_n_steps(
         Stepper stepper , System system , const State &start_state ,
         Time start_time , Time dt , size_t num_of_steps ,
-        Observer observer )
+        Observer observer , StepOverflowChecker checker )
+{
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+    typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+    typedef typename stepper_type::stepper_category stepper_category;
+
+    return detail::integrate_n_steps(
+            checked_stepper<stepper_type, checker_type>(stepper, checker),
+            system , start_state ,
+            start_time , dt , num_of_steps ,
+            checked_observer<observer_type, checker_type>(observer, checker),
+            stepper_category() );
+}
+
+
+/**
+* \brief The same function as above, but without checker.
+*/
+template< class Stepper , class System , class State , class Time , class Observer >
+Time integrate_n_steps(
+        Stepper stepper , System system , State &start_state ,
+        Time start_time , Time dt , size_t num_of_steps , Observer observer )
 {
-    typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
+
     return detail::integrate_n_steps(
-                 stepper , system , start_state ,
-                 start_time , dt , num_of_steps ,
-                 observer , stepper_category() );
+            stepper , system , start_state ,
+            start_time , dt , num_of_steps ,
+            observer , stepper_category() );
 }
 
+/**
+* \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+*/
+template< class Stepper , class System , class State , class Time , class Observer >
+Time integrate_n_steps(
+        Stepper stepper , System system , const State &start_state ,
+        Time start_time , Time dt , size_t num_of_steps , Observer observer )
+{
+    typedef typename odeint::unwrap_reference<Stepper>::type::stepper_category stepper_category;
+
+    return detail::integrate_n_steps(
+            stepper , system , start_state ,
+            start_time , dt , num_of_steps ,
+            observer , stepper_category() );
+}
 
 /**
  * \brief The same function as above, but without observer calls.
@@ -72,7 +117,8 @@ Time integrate_n_steps(
         Stepper stepper , System system , State &start_state ,
         Time start_time , Time dt , size_t num_of_steps )
 {
-    return integrate_n_steps( stepper , system , start_state , start_time , dt , num_of_steps , null_observer() );
+    return integrate_n_steps(stepper, system, start_state, start_time,
+                             dt, num_of_steps, null_observer());
 }
 
 /**
@@ -83,7 +129,8 @@ Time integrate_n_steps(
         Stepper stepper , System system , const State &start_state ,
         Time start_time , Time dt , size_t num_of_steps )
 {
-    return integrate_n_steps( stepper , system , start_state , start_time , dt , num_of_steps , null_observer() );
+    return integrate_n_steps(stepper, system, start_state, start_time,
+                             dt, num_of_steps, null_observer());
 }
 
 
@@ -102,6 +149,11 @@ Time integrate_n_steps(
      * t0 + n*dt. If a DenseOutputStepper is used, the step size also may vary
      * and the dense output is used to call the observer at equidistant time
      * points. The final integration time is always t0 + num_of_steps*dt.
+     * If a max_step_checker is provided as StepOverflowChecker, a
+     * no_progress_errror is thrown if too many steps (default: 500) are
+     * performed without progress, i.e. in between observer calls. If no
+     * checker is provided, no such overflow check is performed.
+
      *
      * \param stepper The stepper to be used for numerical integration.
      * \param system Function/Functor defining the rhs of the ODE.
@@ -111,6 +163,8 @@ Time integrate_n_steps(
      * time step of the integration.
      * \param num_of_steps Number of steps to be performed
      * \param observer Function/Functor called at equidistant time intervals.
+     * \param checker [optional] Functor to check for step count overflows, if no
+     * checker is provided, no exception is thrown.
      * \return The number of steps performed.
      */
 
diff --git a/boost/numeric/odeint/integrate/integrate_times.hpp b/boost/numeric/odeint/integrate/integrate_times.hpp
index c0ecda6c6d..79fba4f1b3 100644
--- a/boost/numeric/odeint/integrate/integrate_times.hpp
+++ b/boost/numeric/odeint/integrate/integrate_times.hpp
@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
 
  Distributed under the Boost Software License, Version 1.0.
  (See accompanying file LICENSE_1_0.txt or
@@ -24,6 +24,7 @@
 
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 #include <boost/numeric/odeint/integrate/null_observer.hpp>
+#include <boost/numeric/odeint/integrate/check_adapter.hpp>
 #include <boost/numeric/odeint/integrate/detail/integrate_times.hpp>
 
 namespace boost {
@@ -32,8 +33,91 @@ namespace odeint {
 
 
 /*
+ * \brief Integrates while calling the observer at the time points given by sequence [times_start, time_end)
  * the two overloads are needed in order to solve the forwarding problem
  */
+template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer , class StepOverflowChecker >
+size_t integrate_times(
+        Stepper stepper , System system , State &start_state ,
+        TimeIterator times_start , TimeIterator times_end , Time dt ,
+        Observer observer , StepOverflowChecker checker )
+{
+    // unwrap references
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+    typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+    typedef typename stepper_type::stepper_category stepper_category;
+
+    // pass on checked stepper and observer
+    // checked_stepper/observer use references internally, so passing by value is fine
+    return detail::integrate_times(
+            checked_stepper<stepper_type, checker_type>(stepper, checker) ,
+            system , start_state ,
+            times_start , times_end , dt ,
+            checked_observer<observer_type, checker_type>(observer, checker),
+            stepper_category() );
+}
+
+/**
+ * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+ */
+template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer , class StepOverflowChecker >
+size_t integrate_times(
+        Stepper stepper , System system , const State &start_state ,
+        TimeIterator times_start , TimeIterator times_end , Time dt ,
+        Observer observer , StepOverflowChecker checker )
+{
+    typedef typename odeint::unwrap_reference< Stepper >::type stepper_type;
+    typedef typename odeint::unwrap_reference< Observer >::type observer_type;
+    typedef typename odeint::unwrap_reference< StepOverflowChecker >::type checker_type;
+    typedef typename stepper_type::stepper_category stepper_category;
+
+    stepper_type &st = stepper;
+    observer_type &obs = observer;
+    checker_type &chk = checker;
+
+    return detail::integrate_times(
+            checked_stepper<stepper_type, checker_type>(stepper, checker) ,
+            system , start_state ,
+            times_start , times_end , dt ,
+            checked_observer<observer_type, checker_type>(observer, checker),
+            stepper_category() );
+}
+
+/**
+ * \brief The same function as above, but with the observation times given as range.
+ */
+template< class Stepper , class System , class State , class TimeRange , class Time , class Observer , class StepOverflowChecker >
+size_t integrate_times(
+        Stepper stepper , System system , State &start_state ,
+        const TimeRange &times , Time dt ,
+        Observer observer , StepOverflowChecker checker )
+{
+    return integrate_times(
+            stepper , system , start_state ,
+            boost::begin( times ) , boost::end( times ) , dt , observer , checker );
+}
+
+/**
+ * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+ */
+template< class Stepper , class System , class State , class TimeRange , class Time , class Observer , class StepOverflowChecker >
+size_t integrate_times(
+        Stepper stepper , System system , const State &start_state ,
+        const TimeRange &times , Time dt ,
+        Observer observer , StepOverflowChecker checker )
+{
+    return integrate_times(
+            stepper , system , start_state ,
+            boost::begin( times ) , boost::end( times ) , dt , observer , checker );
+}
+
+
+
+
+/*
+* The same functions as above, but without a StepOverflowChecker
+*/
 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
 size_t integrate_times(
         Stepper stepper , System system , State &start_state ,
@@ -41,6 +125,7 @@ size_t integrate_times(
         Observer observer )
 {
     typedef typename odeint::unwrap_reference< Stepper >::type::stepper_category stepper_category;
+    // simply don't use checked_* adapters
     return detail::integrate_times(
             stepper , system , start_state ,
             times_start , times_end , dt ,
@@ -48,8 +133,8 @@ size_t integrate_times(
 }
 
 /**
- * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
- */
+* \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+*/
 template< class Stepper , class System , class State , class TimeIterator , class Time , class Observer >
 size_t integrate_times(
         Stepper stepper , System system , const State &start_state ,
@@ -64,8 +149,8 @@ size_t integrate_times(
 }
 
 /**
- * \brief The same function as above, but without observer calls.
- */
+* \brief The same function as above, but with the observation times given as range.
+*/
 template< class Stepper , class System , class State , class TimeRange , class Time , class Observer >
 size_t integrate_times(
         Stepper stepper , System system , State &start_state ,
@@ -78,8 +163,8 @@ size_t integrate_times(
 }
 
 /**
- * \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
- */
+* \brief Solves the forwarding problem, can be called with Boost.Range as start_state.
+*/
 template< class Stepper , class System , class State , class TimeRange , class Time , class Observer >
 size_t integrate_times(
         Stepper stepper , System system , const State &start_state ,
@@ -88,12 +173,10 @@ size_t integrate_times(
 {
     return integrate_times(
             stepper , system , start_state ,
-            boost::begin( times ) , boost::end( times ) , dt , observer );
+            boost::begin( times ) , boost::end( times ) , dt , observer);
 }
 
 
-
-
 /********* DOXYGEN ***********/
 
     /**
@@ -110,6 +193,10 @@ size_t integrate_times(
      * If a DenseOutputStepper is provided, the dense output functionality is
      * used to call the observer at the given times. The end time of the 
      * integration is always *(end_time-1).
+     * If a max_step_checker is provided as StepOverflowChecker, a
+     * no_progress_error is thrown if too many steps (default: 500) are
+     * performed without progress, i.e. in between observer calls. If no
+     * checker is provided, no such overflow check is performed.
      *
      * \param stepper The stepper to be used for numerical integration.
      * \param system Function/Functor defining the rhs of the ODE.
@@ -119,6 +206,8 @@ size_t integrate_times(
      * \param dt The time step between observer calls, _not_ necessarily the 
      * time step of the integration.
      * \param observer Function/Functor called at equidistant time intervals.
+     * \param checker [optional] Functor to check for step count overflows, if no
+     * checker is provided, no exception is thrown.
      * \return The number of steps performed.
      */
 
diff --git a/boost/numeric/odeint/integrate/max_step_checker.hpp b/boost/numeric/odeint/integrate/max_step_checker.hpp
new file mode 100644
index 0000000000..6808a57bdb
--- /dev/null
+++ b/boost/numeric/odeint/integrate/max_step_checker.hpp
@@ -0,0 +1,114 @@
+/*
+ [auto_generated]
+ boost/numeric/odeint/integrate/max_step_checker.hpp
+
+ [begin_description]
+ Throws exception if too many steps are performed.
+ [end_description]
+
+ Copyright 2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+
+
+#ifndef BOOST_NUMERIC_ODEINT_INTEGRATE_MAX_STEP_CHECKER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_INTEGRATE_MAX_STEP_CHECKER_HPP_INCLUDED
+
+#include <stdexcept>
+#include <cstdio>
+
+#include <boost/throw_exception.hpp>
+#include <boost/numeric/odeint/util/odeint_error.hpp>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+/**
+ * \brief A class for performing overflow checks on the step count in integrate functions.
+ *
+ * Provide an instance of this class to integrate functions if you want to throw a runtime error if
+ * too many steps are performed without progress during the integrate routine.
+ */
+class max_step_checker
+{
+public:
+
+protected:
+    const int m_max_steps;
+    int m_steps;
+
+public:
+    /**
+     * \brief Construct the max_step_checker.
+     * max_steps is the maximal number of iterations allowed before runtime_error is thrown.
+     */
+    max_step_checker(const int max_steps = 500)
+        : m_max_steps(max_steps)
+    {
+        reset();
+    }
+
+    /**
+     * \brief Resets the max_step_checker by setting the internal counter to 0.
+     */
+    void reset()
+    {
+        m_steps = 0;
+    }
+
+    /**
+     * \brief Increases the counter and performs the iteration check
+     */
+    void operator()(void)
+    {
+        if( m_steps++ >= m_max_steps )
+        {
+            char error_msg[200];
+            sprintf(error_msg, "Max number of iterations exceeded (%d).", m_max_steps);
+            BOOST_THROW_EXCEPTION( no_progress_error(error_msg) );
+        }
+    }
+};
+
+
+/**
+ * \brief A class for performing overflow checks on the failed step count in step size adjustments.
+ *
+ * Used internally within the dense output stepper and integrate routines.
+ */
+class failed_step_checker : public max_step_checker
+{
+
+public:
+    /**
+     * \brief Construct the failed_step_checker.
+     * max_steps is the maximal number of iterations allowed before runtime_error is thrown.
+     */
+    failed_step_checker(const int max_steps = 500)
+            : max_step_checker(max_steps)
+    {}
+
+    /**
+     * \brief Increases the counter and performs the iteration check
+     */
+    void operator()(void)
+    {
+        if( m_steps++ >= m_max_steps )
+        {
+            char error_msg[200];
+            sprintf(error_msg, "Max number of iterations exceeded (%d). A new step size was not found.", m_max_steps);
+            BOOST_THROW_EXCEPTION( step_adjustment_error(error_msg) );
+        }
+    }
+};
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+#endif
\ No newline at end of file
diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
index e71008a993..4b908333ba 100644
--- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
+++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
@@ -90,9 +90,11 @@ public:
 
     bulirsch_stoer(
         value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
-        value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 )
+        value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
+        time_type max_dt = static_cast<time_type>(0))
         : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() ,
           m_last_step_rejected( false ) , m_first( true ) ,
+          m_max_dt(max_dt) ,
           m_interval_sequence( m_k_max+1 ) ,
           m_coeff( m_k_max+1 ) ,
           m_cost( m_k_max+1 ) ,
@@ -189,6 +191,14 @@ public:
     template< class System , class StateIn , class DerivIn , class StateOut >
     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
     {
+        if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
+        {
+            // given step size is bigger then max_dt
+            // set limit and return fail
+            dt = m_max_dt;
+            return fail;
+        }
+
         BOOST_USING_STD_MIN();
         BOOST_USING_STD_MAX();
 
@@ -311,6 +321,11 @@ public:
 
         if( !m_last_step_rejected || boost::numeric::odeint::detail::less_with_sign(new_h, dt, dt) )
         {
+            // limit step size
+            if( m_max_dt != static_cast<time_type>(0) )
+            {
+                new_h = detail::min_abs(m_max_dt, new_h);
+            }
             m_dt_last = new_h;
             dt = new_h;
         }
@@ -474,6 +489,7 @@ private:
 
     time_type m_dt_last;
     time_type m_t_last;
+    time_type m_max_dt;
 
     size_t m_current_k_opt;
 
@@ -493,7 +509,7 @@ private:
 
     state_table_type m_table; // sequence of states for extrapolation
 
-    const value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
+    value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
 };
 
 
diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
index 3018524ec9..d876ca3d36 100644
--- a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
+++ b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
@@ -6,7 +6,7 @@
  Implementaiton of the Burlish-Stoer method with dense output
  [end_description]
 
- Copyright 2011-2013 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2011-2013 Karsten Ahnert
  Copyright 2012 Christoph Koke
 
@@ -43,6 +43,8 @@
 #include <boost/numeric/odeint/util/resizer.hpp>
 #include <boost/numeric/odeint/util/unit_helper.hpp>
 
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
+
 #include <boost/type_traits.hpp>
 
 
@@ -97,8 +99,10 @@ public:
     bulirsch_stoer_dense_out(
         value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 ,
         value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 ,
+        time_type max_dt = static_cast<time_type>(0) ,
         bool control_interpolation = false )
-        : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , 
+        : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) ,
+          m_max_dt(max_dt) ,
           m_control_interpolation( control_interpolation) ,
           m_last_step_rejected( false ) , m_first( true ) ,
           m_current_state_x1( true ) ,
@@ -149,6 +153,14 @@ public:
     template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , DerivOut &dxdt_new , time_type &dt )
     {
+        if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
+        {
+            // given step size is bigger then max_dt
+            // set limit and return fail
+            dt = m_max_dt;
+            return fail;
+        }
+
         BOOST_USING_STD_MIN();
         BOOST_USING_STD_MAX();
         using std::pow;
@@ -275,7 +287,14 @@ public:
         }
         //set next stepsize
         if( !m_last_step_rejected || (new_h < dt) )
+        {
+            // limit step size
+            if( m_max_dt != static_cast<time_type>(0) )
+            {
+                new_h = detail::min_abs(m_max_dt, new_h);
+            }
             dt = new_h;
+        }
 
         m_last_step_rejected = reject;
         if( reject )
@@ -301,23 +320,20 @@ public:
     template< class System >
     std::pair< time_type , time_type > do_step( System system )
     {
-        const size_t max_count = 1000;
-
         if( m_first )
         {
             typename odeint::unwrap_reference< System >::type &sys = system;
             sys( get_current_state() , get_current_deriv() , m_t );
         }
 
+        failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
         controlled_step_result res = fail;
         m_t_last = m_t;
-        size_t count = 0;
         while( res == fail )
         {
             res = try_step( system , get_current_state() , get_current_deriv() , m_t , get_old_state() , get_old_deriv() , m_dt );
             m_first = false;
-            if( count++ == max_count )
-                throw std::overflow_error( "bulirsch_stoer : too much iterations!");
+            fail_checker();  // check for overflow of failed steps
         }
         toggle_current_state();
         return std::make_pair( m_t_last , m_t );
@@ -412,15 +428,15 @@ private:
         BOOST_USING_STD_MAX();
         using std::pow;
 
-        value_type expo=1.0/(m_interval_sequence[k-1]);
+        value_type expo = static_cast<value_type>(1)/(m_interval_sequence[k-1]);
         value_type facmin = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , expo );
         value_type fac;
         if (error == 0.0)
-            fac=1.0/facmin;
+            fac = static_cast<value_type>(1)/facmin;
         else
         {
             fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo );
-            fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( facmin/STEPFAC4 , min BOOST_PREVENT_MACRO_SUBSTITUTION( 1.0/facmin , fac ) );
+            fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>( facmin/STEPFAC4 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast<value_type>(static_cast<value_type>(1)/facmin) , fac ) );
         }
         return h*fac;
     }
@@ -646,6 +662,8 @@ private:
     default_error_checker< value_type, algebra_type , operations_type > m_error_checker;
     modified_midpoint_dense_out< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint;
 
+    time_type m_max_dt;
+
     bool m_control_interpolation;
 
     bool m_last_step_rejected;
@@ -684,7 +702,7 @@ private:
 
     //wrapped_state_type m_a1 , m_a2 , m_a3 , m_a4;
 
-    const value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
+    value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2;
 };
 
 
diff --git a/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp b/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
index 509192482c..8ae627fe1c 100644
--- a/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
+++ b/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp
@@ -6,7 +6,7 @@
  [end_description]
 
  Copyright 2010-2013 Karsten Ahnert
- Copyright 2010-2013 Mario Mulansky
+ Copyright 2010-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
 
  Distributed under the Boost Software License, Version 1.0.
@@ -33,6 +33,7 @@
 #include <boost/numeric/odeint/util/state_wrapper.hpp>
 #include <boost/numeric/odeint/util/is_resizeable.hpp>
 #include <boost/numeric/odeint/util/resizer.hpp>
+#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
 
 #include <boost/numeric/odeint/algebra/range_algebra.hpp>
 #include <boost/numeric/odeint/algebra/default_operations.hpp>
@@ -64,23 +65,24 @@ public:
             value_type eps_abs = static_cast< value_type >( 1.0e-6 ) ,
             value_type eps_rel = static_cast< value_type >( 1.0e-6 ) ,
             value_type a_x = static_cast< value_type >( 1 ) ,
-            value_type a_dxdt = static_cast< value_type >( 1 ) )
-    : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
+            value_type a_dxdt = static_cast< value_type >( 1 ))
+        : m_eps_abs( eps_abs ) , m_eps_rel( eps_rel ) , m_a_x( a_x ) , m_a_dxdt( a_dxdt )
     { }
 
 
-    template< class State , class Deriv , class Err , class Time >
+    template< class State , class Deriv , class Err, class Time >
     value_type error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
     {
         return error( algebra_type() , x_old , dxdt_old , x_err , dt );
     }
 
-    template< class State , class Deriv , class Err , class Time >
+    template< class State , class Deriv , class Err, class Time >
     value_type error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
     {
+        using std::abs;
         // this overwrites x_err !
         algebra.for_each3( x_err , x_old , dxdt_old ,
-                typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * get_unit_value( dt ) ) );
+                typename operations_type::template rel_error< value_type >( m_eps_abs , m_eps_rel , m_a_x , m_a_dxdt * abs(get_unit_value( dt )) ) );
 
         // value_type res = algebra.reduce( x_err ,
         //        typename operations_type::template maximum< value_type >() , static_cast< value_type >( 0 ) );
@@ -97,9 +99,73 @@ private:
 };
 
 
+template< typename Value, typename Time >
+class default_step_adjuster
+{
+public:
+    typedef Time time_type;
+    typedef Value value_type;
+
+    default_step_adjuster(const time_type max_dt=static_cast<time_type>(0))
+            : m_max_dt(max_dt)
+    {}
 
 
+    time_type decrease_step(time_type dt, const value_type error, const int error_order) const
+    {
+        // returns the decreased time step
+        BOOST_USING_STD_MIN();
+        BOOST_USING_STD_MAX();
+        using std::pow;
 
+        dt *= max
+        BOOST_PREVENT_MACRO_SUBSTITUTION(
+                static_cast<value_type>( static_cast<value_type>(9) / static_cast<value_type>(10) *
+                                         pow(error, static_cast<value_type>(-1) / (error_order - 1))),
+                static_cast<value_type>( static_cast<value_type>(1) / static_cast<value_type> (5)));
+        if(m_max_dt != static_cast<time_type >(0))
+            // limit to maximal stepsize even when decreasing
+            dt = detail::min_abs(dt, m_max_dt);
+        return dt;
+    }
+
+    time_type increase_step(time_type dt, value_type error, const int stepper_order) const
+    {
+        // returns the increased time step
+        BOOST_USING_STD_MIN();
+        BOOST_USING_STD_MAX();
+        using std::pow;
+
+        // adjust the size if dt is smaller than max_dt (providede max_dt is not zero)
+        if(error < 0.5)
+        {
+            // error should be > 0
+            error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
+                    static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) ,
+                    error);
+            time_type dt_old = dt;
+            //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
+            dt *= static_cast<value_type>(9)/static_cast<value_type>(10) *
+                  pow(error, static_cast<value_type>(-1) / stepper_order);
+            if(m_max_dt != static_cast<time_type >(0))
+                // limit to maximal stepsize
+                dt = detail::min_abs(dt, m_max_dt);
+        }
+        return dt;
+    }
+
+    bool check_step_size_limit(const time_type dt)
+    {
+        if(m_max_dt != static_cast<time_type >(0))
+            return detail::less_eq_with_sign(dt, m_max_dt, dt);
+        return true;
+    }
+
+    time_type get_max_dt() { return m_max_dt; }
+
+private:
+    time_type m_max_dt;
+};
 
 
 
@@ -109,8 +175,10 @@ private:
 template<
 class ErrorStepper ,
 class ErrorChecker = default_error_checker< typename ErrorStepper::value_type ,
-typename ErrorStepper::algebra_type ,
-typename ErrorStepper::operations_type > ,
+    typename ErrorStepper::algebra_type ,
+    typename ErrorStepper::operations_type > ,
+class StepAdjuster = default_step_adjuster< typename ErrorStepper::value_type ,
+    typename ErrorStepper::time_type > ,
 class Resizer = typename ErrorStepper::resizer_type ,
 class ErrorStepperCategory = typename ErrorStepper::stepper_category
 >
@@ -139,11 +207,13 @@ class controlled_runge_kutta ;
  * \tparam Resizer The resizer policy type.
  */
 template<
-class ErrorStepper ,
-class ErrorChecker ,
+class ErrorStepper,
+class ErrorChecker,
+class StepAdjuster,
 class Resizer
 >
-class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag >
+class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster, Resizer ,
+        explicit_error_stepper_tag >
 {
 
 public:
@@ -157,13 +227,15 @@ public:
     typedef typename stepper_type::operations_type operations_type;
     typedef Resizer resizer_type;
     typedef ErrorChecker error_checker_type;
+    typedef StepAdjuster step_adjuster_type;
     typedef explicit_controlled_stepper_tag stepper_category;
 
 #ifndef DOXYGEN_SKIP
     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
 
-    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
+    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster ,
+            Resizer , explicit_error_stepper_tag > controlled_stepper_type;
 #endif //DOXYGEN_SKIP
 
 
@@ -174,9 +246,10 @@ public:
      */
     controlled_runge_kutta(
             const error_checker_type &error_checker = error_checker_type( ) ,
+            const step_adjuster_type &step_adjuster = step_adjuster_type() ,
             const stepper_type &stepper = stepper_type( )
     )
-    : m_stepper( stepper ) , m_error_checker( error_checker )
+        : m_stepper(stepper), m_error_checker(error_checker) , m_step_adjuster(step_adjuster)
     { }
 
 
@@ -338,59 +411,35 @@ public:
     template< class System , class StateIn , class DerivIn , class StateOut >
     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt )
     {
-        BOOST_USING_STD_MIN();
-        BOOST_USING_STD_MAX();
-        using std::pow;
+        if( !m_step_adjuster.check_step_size_limit(dt) )
+        {
+            // given dt was above step size limit - adjust and return fail;
+            dt = m_step_adjuster.get_max_dt();
+            return fail;
+        }
 
         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
 
         // do one step with error calculation
         m_stepper.do_step( system , in , dxdt , t , out , dt , m_xerr.m_v );
 
-        m_max_rel_error = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt );
+        value_type max_rel_err = m_error_checker.error( m_stepper.algebra() , in , dxdt , m_xerr.m_v , dt );
 
-        if( m_max_rel_error > 1.0 )
+        if( max_rel_err > 1.0 )
         {
-            // error too large - decrease dt ,limit scaling factor to 0.2 and reset state
-            dt *= max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) *
-                                                         pow( m_max_rel_error , static_cast<value_type>(-1) / ( m_stepper.error_order() - 1 ) ) ) ,
-                                                         static_cast<value_type>( static_cast<value_type>(1)/static_cast<value_type> (5) ) );
+            // error too big, decrease step size and reject this step
+            dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
             return fail;
-        }
-        else
+        } else
         {
-            if( m_max_rel_error < 0.5 )
-            {
-                // error should be > 0
-                m_max_rel_error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
-                            static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(m_stepper.stepper_order()) ) ) ,
-                            m_max_rel_error );
-                //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
-                t += dt;
-                dt *= static_cast<value_type>(9)/static_cast<value_type>(10) * pow( m_max_rel_error ,
-                                                                                    static_cast<value_type>(-1) / m_stepper.stepper_order() );
-                return success;
-            }
-            else
-            {
-                t += dt;
-                return success;
-            }
+            // otherwise, increase step size and accept
+            t += dt;
+            dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
+            return success;
         }
     }
 
     /**
-     * \brief Returns the error of the last step.
-     *
-     * returns The last error of the step.
-     */
-    value_type last_error( void ) const
-    {
-        return m_max_rel_error;
-    }
-
-
-    /**
      * \brief Adjust the size of all temporaries in the stepper manually.
      * \param x A state from which the size of the temporaries to be resized is deduced.
      */
@@ -455,6 +504,7 @@ private:
 
     stepper_type m_stepper;
     error_checker_type m_error_checker;
+    step_adjuster_type m_step_adjuster;
 
     resizer_type m_dxdt_resizer;
     resizer_type m_xerr_resizer;
@@ -463,7 +513,6 @@ private:
     wrapped_deriv_type m_dxdt;
     wrapped_state_type m_xerr;
     wrapped_state_type m_xnew;
-    value_type m_max_rel_error;
 };
 
 
@@ -498,9 +547,10 @@ private:
 template<
 class ErrorStepper ,
 class ErrorChecker ,
+class StepAdjuster ,
 class Resizer
 >
-class controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_fsal_tag >
+class controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_fsal_tag >
 {
 
 public:
@@ -514,13 +564,14 @@ public:
     typedef typename stepper_type::operations_type operations_type;
     typedef Resizer resizer_type;
     typedef ErrorChecker error_checker_type;
+    typedef StepAdjuster step_adjuster_type;
     typedef explicit_controlled_stepper_fsal_tag stepper_category;
 
 #ifndef DOXYGEN_SKIP
     typedef typename stepper_type::wrapped_state_type wrapped_state_type;
     typedef typename stepper_type::wrapped_deriv_type wrapped_deriv_type;
 
-    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
+    typedef controlled_runge_kutta< ErrorStepper , ErrorChecker , StepAdjuster , Resizer , explicit_error_stepper_tag > controlled_stepper_type;
 #endif // DOXYGEN_SKIP
 
     /**
@@ -530,9 +581,10 @@ public:
      */
     controlled_runge_kutta(
             const error_checker_type &error_checker = error_checker_type() ,
+            const step_adjuster_type &step_adjuster = step_adjuster_type() ,
             const stepper_type &stepper = stepper_type()
     )
-    : m_stepper( stepper ) , m_error_checker( error_checker ) ,
+    : m_stepper( stepper ) , m_error_checker( error_checker ) , m_step_adjuster(step_adjuster) ,
       m_first_call( true )
     { }
 
@@ -699,10 +751,12 @@ public:
     controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type &t ,
             StateOut &out , DerivOut &dxdt_out , time_type &dt )
     {
-        BOOST_USING_STD_MIN();
-        BOOST_USING_STD_MAX();
-
-        using std::pow;
+        if( !m_step_adjuster.check_step_size_limit(dt) )
+        {
+            // given dt was above step size limit - adjust and return fail;
+            dt = m_step_adjuster.get_max_dt();
+            return fail;
+        }
 
         m_xerr_resizer.adjust_size( in , detail::bind( &controlled_runge_kutta::template resize_m_xerr_impl< StateIn > , detail::ref( *this ) , detail::_1 ) );
 
@@ -715,29 +769,14 @@ public:
 
         if( max_rel_err > 1.0 )
         {
-            // error too large - decrease dt ,limit scaling factor to 0.2 and reset state
-            dt *= max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) *
-                                                                 pow( max_rel_err , static_cast<value_type>(-1) / ( m_stepper.error_order() - 1 ) ) ) ,
-                                                         static_cast<value_type>( static_cast<value_type>(1)/static_cast<value_type> (5)) );
+            // error too big, decrease step size and reject this step
+            dt = m_step_adjuster.decrease_step(dt, max_rel_err, m_stepper.error_order());
             return fail;
         }
-        else
-        {
-            if( max_rel_err < 0.5 )
-            {                //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
-                // error should be > 0
-                max_rel_err = max BOOST_PREVENT_MACRO_SUBSTITUTION ( static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(m_stepper.stepper_order()) ) ) ,
-                                                                     max_rel_err );
-                t += dt;
-                dt *= static_cast<value_type>( static_cast<value_type>(9)/static_cast<value_type>(10) * pow( max_rel_err , static_cast<value_type>(-1) / m_stepper.stepper_order() ) );
-                return success;
-            }
-            else
-            {
-                t += dt;
-                return success;
-            }
-        }
+        // otherwise, increase step size and accept
+        t += dt;
+        dt = m_step_adjuster.increase_step(dt, max_rel_err, m_stepper.stepper_order());
+        return success;
     }
 
 
@@ -863,6 +902,7 @@ private:
 
     stepper_type m_stepper;
     error_checker_type m_error_checker;
+    step_adjuster_type m_step_adjuster;
 
     resizer_type m_dxdt_resizer;
     resizer_type m_xerr_resizer;
@@ -890,12 +930,14 @@ private:
  * It is used by the controlled_runge_kutta steppers.
  *
  * \tparam Value The value type.
+ * \tparam Time The time type.
  * \tparam Algebra The algebra type.
  * \tparam Operations The operations type.
  */
 
     /**
-     * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt )
+     * \fn default_error_checker( value_type eps_abs , value_type eps_rel , value_type a_x , value_type a_dxdt ,
+     * time_type max_dt)
      * \brief Constructs the error checker.
      *
      * The error is calculated as follows: ???? 
@@ -904,10 +946,11 @@ private:
      * \param eps_rel Relative tolerance level.
      * \param a_x Factor for the weight of the state.
      * \param a_dxdt Factor for the weight of the derivative.
+     * \param max_dt Maximum allowed step size.
      */
     
     /**
-     * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
+     * \fn error( const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
      * \brief Calculates the error level.
      *
      * If the returned error level is greater than 1, the estimated error was
@@ -922,7 +965,7 @@ private:
      */
 
     /**
-     * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , Time dt ) const
+     * \fn error( algebra_type &algebra , const State &x_old , const Deriv &dxdt_old , Err &x_err , time_type dt ) const
      * \brief Calculates the error level using a given algebra.
      *
      * If the returned error level is greater than 1, the estimated error was
@@ -937,6 +980,31 @@ private:
      * \return error
      */
 
+    /**
+     * \fn time_type decrease_step(const time_type dt, const value_type error, const int error_order)
+     * \brief Returns a decreased step size based on the given error and order
+     *
+     * Calculates a new smaller step size based on the given error and its order.
+     *
+     * \param dt The old step size.
+     * \param error The computed error estimate.
+     * \param error_order The error order of the stepper.
+     * \return dt_new The new, reduced step size.
+     */
+
+    /**
+     * \fn time_type increase_step(const time_type dt, const value_type error, const int error_order)
+     * \brief Returns an increased step size based on the given error and order.
+     *
+     * Calculates a new bigger step size based on the given error and its order. If max_dt != 0, the
+     * new step size is limited to max_dt.
+     *
+     * \param dt The old step size.
+     * \param error The computed error estimate.
+     * \param error_order The order of the stepper.
+     * \return dt_new The new, increased step size.
+     */
+
 
 } // odeint
 } // numeric
diff --git a/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp b/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp
index b685d92bcf..94abc5af99 100644
--- a/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp
+++ b/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp
@@ -8,7 +8,7 @@
  [end_description]
 
  Copyright 2011-2013 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
 
  Distributed under the Boost Software License, Version 1.0.
@@ -37,6 +37,8 @@
 #include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
+
 namespace boost {
 namespace numeric {
 namespace odeint {
@@ -206,6 +208,15 @@ public:
         return m_t_old;
     }
 
+    /**
+     * \brief Returns the current time step.
+     * \return dt.
+     */
+    time_type current_time_step( void ) const
+    {
+        return m_dt;
+    }
+
 
 private:
 
@@ -312,8 +323,6 @@ public:
     template< class System >
     std::pair< time_type , time_type > do_step( System system )
     {
-        const size_t max_count = 1000;
-
         if( !m_is_deriv_initialized )
         {
             typename odeint::unwrap_reference< System >::type &sys = system;
@@ -321,15 +330,14 @@ public:
             m_is_deriv_initialized = true;
         }
 
+        failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
         controlled_step_result res = fail;
         m_t_old = m_t;
-        size_t count = 0;
         do
         {
             res = m_stepper.try_step( system , get_current_state() , get_current_deriv() , m_t ,
                                       get_old_state() , get_old_deriv() , m_dt );
-            if( count++ == max_count )
-                BOOST_THROW_EXCEPTION( std::overflow_error( "dense_output_controlled_explicit : too much iterations!") );
+            fail_checker();  // check for overflow of failed steps
         }
         while( res == fail );
         toggle_current_state();
diff --git a/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp b/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp
index 54b7746c6e..bad2bedfef 100644
--- a/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp
+++ b/boost/numeric/odeint/stepper/generation/generation_controlled_runge_kutta.hpp
@@ -34,15 +34,23 @@ struct controller_factory< Stepper , controlled_runge_kutta< Stepper > >
     typedef Stepper stepper_type;
     typedef controlled_runge_kutta< stepper_type > controller_type;
     typedef typename controller_type::error_checker_type error_checker_type;
+    typedef typename controller_type::step_adjuster_type step_adjuster_type;
     typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::value_type time_type;
 
     controller_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper )
     {
-        return controller_type( error_checker_type( abs_error , rel_error ) , stepper );
+        return controller_type( error_checker_type( abs_error , rel_error ) ,
+                                step_adjuster_type() , stepper );
     }
-};
-
 
+    controller_type operator()( value_type abs_error , value_type rel_error ,
+                                time_type max_dt, const stepper_type &stepper )
+    {
+        return controller_type( error_checker_type( abs_error , rel_error ) ,
+                                step_adjuster_type(max_dt) , stepper );
+    }
+};
 
 
 } // odeint
diff --git a/boost/numeric/odeint/stepper/generation/generation_dense_output_runge_kutta.hpp b/boost/numeric/odeint/stepper/generation/generation_dense_output_runge_kutta.hpp
index 214d09c1a4..276358a12e 100644
--- a/boost/numeric/odeint/stepper/generation/generation_dense_output_runge_kutta.hpp
+++ b/boost/numeric/odeint/stepper/generation/generation_dense_output_runge_kutta.hpp
@@ -33,12 +33,23 @@ struct dense_output_factory< Stepper , dense_output_runge_kutta< controlled_rung
     typedef Stepper stepper_type;
     typedef controlled_runge_kutta< stepper_type > controller_type;
     typedef typename controller_type::error_checker_type error_checker_type;
+    typedef typename controller_type::step_adjuster_type step_adjuster_type;
     typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::time_type time_type;
     typedef dense_output_runge_kutta< controller_type > dense_output_type;
 
     dense_output_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper )
     {
-        return dense_output_type( controller_type( error_checker_type( abs_error , rel_error ) , stepper ) );
+        return dense_output_type( controller_type( error_checker_type( abs_error , rel_error ) ,
+                                                   step_adjuster_type() , stepper ) );
+    }
+
+    dense_output_type operator()( value_type abs_error , value_type rel_error ,
+                                  time_type max_dt , const stepper_type &stepper )
+    {
+        return dense_output_type(
+                controller_type( error_checker_type( abs_error , rel_error) ,
+                                 step_adjuster_type( max_dt ) , stepper ) );
     }
 };
 
diff --git a/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp b/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp
index d5322c8d73..366cb5e9c4 100644
--- a/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp
+++ b/boost/numeric/odeint/stepper/generation/generation_rosenbrock4.hpp
@@ -54,12 +54,19 @@ struct dense_output_factory< Stepper , rosenbrock4_dense_output< rosenbrock4_con
     typedef Stepper stepper_type;
     typedef rosenbrock4_controller< stepper_type > controller_type;
     typedef typename stepper_type::value_type value_type;
+    typedef typename stepper_type::time_type time_type;
     typedef rosenbrock4_dense_output< controller_type > dense_output_type;
 
     dense_output_type operator()( value_type abs_error , value_type rel_error , const stepper_type &stepper )
     {
         return dense_output_type( controller_type( abs_error , rel_error , stepper ) );
     }
+
+    dense_output_type operator()( value_type abs_error , value_type rel_error ,
+                                  time_type max_dt, const stepper_type &stepper )
+    {
+        return dense_output_type( controller_type( abs_error , rel_error , max_dt , stepper ) );
+    }
 };
 
 
diff --git a/boost/numeric/odeint/stepper/generation/make_controlled.hpp b/boost/numeric/odeint/stepper/generation/make_controlled.hpp
index 603978fec8..61bc8f1920 100644
--- a/boost/numeric/odeint/stepper/generation/make_controlled.hpp
+++ b/boost/numeric/odeint/stepper/generation/make_controlled.hpp
@@ -43,6 +43,15 @@ struct controller_factory
     {
         return Controller( abs_error , rel_error , stepper );
     }
+
+    Controller operator()(
+            typename Stepper::value_type abs_error ,
+            typename Stepper::value_type rel_error ,
+            typename Stepper::time_type max_dt ,
+            const Stepper &stepper )
+    {
+        return Controller( abs_error , rel_error , max_dt, stepper );
+    }
 };
 
 
@@ -72,6 +81,19 @@ typename result_of::make_controlled< Stepper >::type make_controlled(
 }
 
 
+template< class Stepper >
+typename result_of::make_controlled< Stepper >::type make_controlled(
+        typename Stepper::value_type abs_error ,
+        typename Stepper::value_type rel_error ,
+        typename Stepper::time_type max_dt ,
+        const Stepper & stepper = Stepper() )
+{
+    typedef Stepper stepper_type;
+    typedef typename result_of::make_controlled< stepper_type >::type controller_type;
+    typedef controller_factory< stepper_type , controller_type > factory_type;
+    factory_type factory;
+    return factory( abs_error , rel_error , max_dt, stepper );
+}
 
 } // odeint
 } // numeric
diff --git a/boost/numeric/odeint/stepper/generation/make_dense_output.hpp b/boost/numeric/odeint/stepper/generation/make_dense_output.hpp
index a28e31b5d2..fff3590525 100644
--- a/boost/numeric/odeint/stepper/generation/make_dense_output.hpp
+++ b/boost/numeric/odeint/stepper/generation/make_dense_output.hpp
@@ -39,6 +39,15 @@ struct dense_output_factory
     {
         return DenseOutput( abs_error , rel_error , stepper );
     }
+
+    DenseOutput operator()(
+            typename Stepper::value_type abs_error ,
+            typename Stepper::value_type rel_error ,
+            typename Stepper::time_type max_dt ,
+            const Stepper &stepper )
+    {
+        return DenseOutput( abs_error , rel_error , max_dt , stepper );
+    }
 };
 
 
@@ -68,6 +77,19 @@ typename result_of::make_dense_output< Stepper >::type make_dense_output(
 }
 
 
+template< class Stepper >
+typename result_of::make_dense_output< Stepper >::type make_dense_output(
+        typename Stepper::value_type abs_error ,
+        typename Stepper::value_type rel_error ,
+        typename Stepper::time_type max_dt ,
+        const Stepper &stepper = Stepper() )
+{
+    typedef Stepper stepper_type;
+    typedef typename result_of::make_dense_output< stepper_type >::type dense_output_type;
+    typedef dense_output_factory< stepper_type , dense_output_type > factory_type;
+    factory_type factory;
+    return factory( abs_error , rel_error , max_dt, stepper );
+}
 
 
 } // odeint
diff --git a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
index df4e6c48b6..0e5edd32c0 100644
--- a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
+++ b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
@@ -27,6 +27,7 @@
 
 #include <boost/numeric/odeint/util/copy.hpp>
 #include <boost/numeric/odeint/util/is_resizeable.hpp>
+#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
 
 #include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
 
@@ -55,12 +56,20 @@ public:
     typedef rosenbrock4_controller< Stepper > controller_type;
 
 
-    rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 , const stepper_type &stepper = stepper_type() )
-    : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) ,
-      m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
-      m_last_rejected( false )
+    rosenbrock4_controller( value_type atol = 1.0e-6 , value_type rtol = 1.0e-6 ,
+                            const stepper_type &stepper = stepper_type() )
+        : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) ,
+          m_max_dt( static_cast<time_type>(0) ) ,
+          m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
+          m_last_rejected( false )
     { }
 
+    rosenbrock4_controller( value_type atol, value_type rtol, time_type max_dt,
+                            const stepper_type &stepper = stepper_type() )
+            : m_stepper( stepper ) , m_atol( atol ) , m_rtol( rtol ) , m_max_dt( max_dt ) ,
+              m_first_step( true ) , m_err_old( 0.0 ) , m_dt_old( 0.0 ) ,
+              m_last_rejected( false )
+    { }
 
     value_type error( const state_type &x , const state_type &xold , const state_type &xerr )
     {
@@ -104,6 +113,14 @@ public:
     boost::numeric::odeint::controlled_step_result
     try_step( System sys , const state_type &x , time_type &t , state_type &xout , time_type &dt )
     {
+        if( m_max_dt != static_cast<time_type>(0) && detail::less_with_sign(m_max_dt, dt, dt) )
+        {
+            // given step size is bigger then max_dt
+            // set limit and return fail
+            dt = m_max_dt;
+            return fail;
+        }
+
         BOOST_USING_STD_MIN();
         BOOST_USING_STD_MAX();
         using std::pow;
@@ -142,7 +159,11 @@ public:
                 min BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) :
                 max BOOST_PREVENT_MACRO_SUBSTITUTION ( dt_new , dt ) );
             t += dt;
-            dt = dt_new;
+            // limit step size to max_dt
+            if( m_max_dt != static_cast<time_type>(0) )
+            {
+                dt = detail::min_abs(m_max_dt, dt_new);
+            }
             m_last_rejected = false;
             return success;
         }
@@ -198,6 +219,7 @@ private:
     wrapped_state_type m_xerr;
     wrapped_state_type m_xnew;
     value_type m_atol , m_rtol;
+    time_type m_max_dt;
     bool m_first_step;
     value_type m_err_old , m_dt_old;
     bool m_last_rejected;
diff --git a/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp b/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp
index 8a866c562e..6695ba6a97 100644
--- a/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp
+++ b/boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp
@@ -7,7 +7,7 @@
  [end_description]
 
  Copyright 2011-2012 Karsten Ahnert
- Copyright 2011-2012 Mario Mulansky
+ Copyright 2011-2015 Mario Mulansky
  Copyright 2012 Christoph Koke
 
  Distributed under the Boost Software License, Version 1.0.
@@ -27,6 +27,8 @@
 #include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>
 #include <boost/numeric/odeint/util/is_resizeable.hpp>
 
+#include <boost/numeric/odeint/integrate/max_step_checker.hpp>
+
 
 namespace boost {
 namespace numeric {
@@ -73,16 +75,13 @@ public:
     template< class System >
     std::pair< time_type , time_type > do_step( System system )
     {
-        const size_t max_count = 1000;
-
+        failed_step_checker fail_checker;  // to throw a runtime_error if step size adjustment fails
         controlled_step_result res = fail;
         m_t_old = m_t;
-        size_t count = 0;
         do
         {
             res = m_stepper.try_step( system , get_current_state() , m_t , get_old_state() , m_dt );
-            if( count++ == max_count )
-                throw std::overflow_error( "rosenbrock4 : too much iterations!");
+            fail_checker();  // check for overflow of failed steps
         }
         while( res == fail );
         m_stepper.stepper().prepare_dense_output();
diff --git a/boost/numeric/odeint/util/detail/less_with_sign.hpp b/boost/numeric/odeint/util/detail/less_with_sign.hpp
index d90e12f7b4..3ffa7ca179 100644
--- a/boost/numeric/odeint/util/detail/less_with_sign.hpp
+++ b/boost/numeric/odeint/util/detail/less_with_sign.hpp
@@ -6,7 +6,7 @@
  Helper function to compare times taking into account the sign of dt
  [end_description]
 
- Copyright 2012-2013 Mario Mulansky
+ Copyright 2012-2015 Mario Mulansky
  Copyright 2012 Karsten Ahnert
 
  Distributed under the Boost Software License, Version 1.0.
@@ -57,7 +57,7 @@ T min_abs( T t1 , T t2 )
 {
     BOOST_USING_STD_MIN();
     BOOST_USING_STD_MAX();
-    if( t1>0 )
+    if( get_unit_value(t1)>0 )
         return min BOOST_PREVENT_MACRO_SUBSTITUTION ( t1 , t2 );
     else
         return max BOOST_PREVENT_MACRO_SUBSTITUTION ( t1 , t2 );
@@ -68,7 +68,7 @@ T max_abs( T t1 , T t2 )
 {
     BOOST_USING_STD_MIN();
     BOOST_USING_STD_MAX();
-    if( t1>0 )
+    if( get_unit_value(t1)>0 )
         return max BOOST_PREVENT_MACRO_SUBSTITUTION ( t1 , t2 );
     else
         return min BOOST_PREVENT_MACRO_SUBSTITUTION ( t1 , t2 );
diff --git a/boost/numeric/odeint/util/odeint_error.hpp b/boost/numeric/odeint/util/odeint_error.hpp
new file mode 100644
index 0000000000..312e5b1afa
--- /dev/null
+++ b/boost/numeric/odeint/util/odeint_error.hpp
@@ -0,0 +1,77 @@
+/*
+ [auto_generated]
+ boost/numeric/odeint/util/odeint_error.hpp
+
+ [begin_description]
+ Runtime Exceptions thrown by odeint
+ [end_description]
+
+ Copyright 2015 Mario Mulansky
+
+ Distributed under the Boost Software License, Version 1.0.
+ (See accompanying file LICENSE_1_0.txt or
+ copy at http://www.boost.org/LICENSE_1_0.txt)
+*/
+
+
+#ifndef BOOST_NUMERIC_ODEINT_UTIL_ODEINT_ERROR_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_UTIL_ODEINT_ERROR_HPP_INCLUDED
+
+#include <stdexcept>
+#include <string>
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+/**
+ * \brief Runtime error thrown by odeint
+ */
+class odeint_error : public std::runtime_error
+{
+public:
+    odeint_error(const std::string &s)
+            : std::runtime_error(s)
+    { }
+};
+
+
+/**
+ * \brief Runtime error thrown from integrate routines
+ *
+ * This Error occures when too many iterations are performed in between two
+ * observer calls in the integrate routines.
+ */
+class no_progress_error : public odeint_error
+{
+public:
+    no_progress_error(const std::string &s)
+            : odeint_error(s)
+    { }
+};
+
+
+/**
+ * \brief Runtime error thrown during stepsize adjustment
+ *
+ * This Error occures when too many iterations are performed without finding
+ * an appropriate new step size. This usually indicates non-continuous points
+ * in the ODE.
+ */
+class step_adjustment_error : public odeint_error
+{
+public:
+    step_adjustment_error(const std::string &s)
+            : odeint_error(s)
+    { }
+};
+
+}
+}
+}
+
+
+
+#endif // BOOST_NUMERIC_ODEINT_UTIL_ODEINT_ERROR_HPP_INCLUDED
diff --git a/boost/numeric/odeint/util/resize.hpp b/boost/numeric/odeint/util/resize.hpp
index 75e9eff448..3645782aee 100644
--- a/boost/numeric/odeint/util/resize.hpp
+++ b/boost/numeric/odeint/util/resize.hpp
@@ -83,7 +83,7 @@ namespace detail {
         }
 
         template< class StateOut , class StateIn >
-        void resize_op( StateOut &x1 , const StateIn &x2 , boost::false_type ) const
+        void resize_op( StateOut &/*x1*/ , const StateIn &/*x2*/ , boost::false_type ) const
         {
         }
 
diff --git a/boost/numeric/odeint/util/resizer.hpp b/boost/numeric/odeint/util/resizer.hpp
index fd969ec3b9..cd27990bda 100644
--- a/boost/numeric/odeint/util/resizer.hpp
+++ b/boost/numeric/odeint/util/resizer.hpp
@@ -79,7 +79,7 @@ struct initially_resizer
 struct never_resizer
 {
     template< class State , class ResizeFunction >
-    bool adjust_size( const State &x , ResizeFunction f )
+    bool adjust_size( const State &/*x*/ , ResizeFunction /*f*/ )
     {
         return false;
     }
diff --git a/boost/numeric/odeint/util/same_size.hpp b/boost/numeric/odeint/util/same_size.hpp
index 8fa6c7a3a0..81ae249d1f 100644
--- a/boost/numeric/odeint/util/same_size.hpp
+++ b/boost/numeric/odeint/util/same_size.hpp
@@ -83,7 +83,7 @@ struct same_size_fusion
     }
 
     template< class S1 , class S2 >
-    bool same_size_op( const S1 &x1 , const S2 &x2 , boost::false_type ) const
+    bool same_size_op( const S1 &/*x1*/ , const S2 &/*x2*/ , boost::false_type ) const
     {
         return true;
     }
diff --git a/boost/numeric/odeint/util/unwrap_reference.hpp b/boost/numeric/odeint/util/unwrap_reference.hpp
index bc7d423a09..e160878391 100644
--- a/boost/numeric/odeint/util/unwrap_reference.hpp
+++ b/boost/numeric/odeint/util/unwrap_reference.hpp
@@ -31,7 +31,7 @@
 namespace boost {
 
 #if BOOST_NUMERIC_ODEINT_CXX11
-template<typename T> struct reference_wrapper;
+template<typename T> class reference_wrapper;
 
 template<typename T> struct unwrap_reference;
 #endif