Imported Upstream version 1.58.0upstream/1.58.0

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

33 files changed, 1015 insertions, 283 deletions

diff --git a/boost/numeric/odeint/algebra/algebra_dispatcher.hpp b/boost/numeric/odeint/algebra/algebra_dispatcher.hpp
index df5aabac96..88cf159e03 100644
--- a/boost/numeric/odeint/algebra/algebra_dispatcher.hpp
+++ b/boost/numeric/odeint/algebra/algebra_dispatcher.hpp
@@ -17,8 +17,9 @@
 #ifndef BOOST_NUMERIC_ODEINT_ALGEBRA_ALGEBRA_DISPATCHER_HPP_INCLUDED
 #define BOOST_NUMERIC_ODEINT_ALGEBRA_ALGEBRA_DISPATCHER_HPP_INCLUDED
 
-#include <complex>
+#include <boost/numeric/odeint/config.hpp>
 
+#include <complex>
 #include <boost/type_traits/is_floating_point.hpp>
 
 #include <boost/numeric/ublas/vector.hpp>
@@ -38,14 +39,14 @@ namespace odeint {
 template< class StateType , class Enabler = void >
 struct algebra_dispatcher_sfinae
 {
-    //  range_algebra is the standard algebra^
+    // range_algebra is the standard algebra
     typedef range_algebra algebra_type;
 };
 
 template< class StateType >
 struct algebra_dispatcher : algebra_dispatcher_sfinae< StateType > { };
 
-//specialize for array
+// specialize for array
 template< class T , size_t N >
 struct algebra_dispatcher< boost::array< T , N > >
 {
@@ -84,4 +85,26 @@ struct algebra_dispatcher< boost::numeric::ublas::matrix< T , L , A > >
 }
 }
 
+#ifdef BOOST_NUMERIC_ODEINT_CXX11
+
+// c++11 mode: specialization for std::array if available
+
+#include <array>
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+    
+// specialize for std::array
+template< class T , size_t N >
+struct algebra_dispatcher< std::array< T , N > >
+{
+    typedef array_algebra algebra_type;
+};
+
+} } }
+
+#endif
+
+
 #endif
diff --git a/boost/numeric/odeint/algebra/array_algebra.hpp b/boost/numeric/odeint/algebra/array_algebra.hpp
index c565d6dd76..471e866fe5 100644
--- a/boost/numeric/odeint/algebra/array_algebra.hpp
+++ b/boost/numeric/odeint/algebra/array_algebra.hpp
@@ -3,7 +3,11 @@
  boost/numeric/odeint/algebra/array_algebra.hpp
 
  [begin_description]
- Algebra for boost::array. Highly specialized for odeint. Const arguments are introduce to work with odeint.
+ Algebra for Arrays. Highly specialized for odeint. Const arguments are
+ introduce to work with odeint.
+ The Array algebra can be used for Array structures with two template
+ parameters: 
+ Array<T, N>
  [end_description]
 
  Copyright 2011-2013 Mario Mulansky
@@ -29,234 +33,251 @@ namespace odeint {
 
 struct array_algebra
 {
-    template< typename T , size_t dim , class Op >
-    static void for_each1( boost::array< T , dim > &s1 , Op op )
+    //template< typename T , size_t dim , class Op >
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each1( Array< T, dim > &s1, Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] );
     }
 
-
-    template< typename T1 , typename T2 , size_t dim , class Op >
-    static void for_each2( boost::array< T1 , dim > &s1 ,
-            const boost::array< T2 , dim > &s2 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each2( Array< T, dim > &s1, const Array< T, dim > &s2,
+                           Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each3( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each3( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] );
     }
 
     /* different const signature - required for the scale_sum_swap2 operation */
-    template< typename T , size_t dim , class Op >
-    static void for_each3( boost::array< T , dim > &s1 ,
-            boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each3( Array< T , dim > &s1 ,
+            Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each4( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each4( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each5( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each5( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each6( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each6( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each7( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each7( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each8( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each8( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each9( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each9( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each10( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each10( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each11( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each11( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each12( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 ,
-            const boost::array< T , dim > &s12 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each12( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 ,
+            const Array< T , dim > &s12 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each13( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 ,
-            const boost::array< T , dim > &s12 ,
-            const boost::array< T , dim > &s13 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each13( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 ,
+            const Array< T , dim > &s12 ,
+            const Array< T , dim > &s13 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each14( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 ,
-            const boost::array< T , dim > &s12 ,
-            const boost::array< T , dim > &s13 ,
-            const boost::array< T , dim > &s14 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each14( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 ,
+            const Array< T , dim > &s12 ,
+            const Array< T , dim > &s13 ,
+            const Array< T , dim > &s14 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] , s14[i] );
     }
 
-    template< typename T , size_t dim , class Op >
-    static void for_each15( boost::array< T , dim > &s1 ,
-            const boost::array< T , dim > &s2 ,
-            const boost::array< T , dim > &s3 ,
-            const boost::array< T , dim > &s4 ,
-            const boost::array< T , dim > &s5 ,
-            const boost::array< T , dim > &s6 ,
-            const boost::array< T , dim > &s7 ,
-            const boost::array< T , dim > &s8 ,
-            const boost::array< T , dim > &s9 ,
-            const boost::array< T , dim > &s10 ,
-            const boost::array< T , dim > &s11 ,
-            const boost::array< T , dim > &s12 ,
-            const boost::array< T , dim > &s13 ,
-            const boost::array< T , dim > &s14 ,
-            const boost::array< T , dim > &s15 , Op op )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim, class Op >
+    static void for_each15( Array< T , dim > &s1 ,
+            const Array< T , dim > &s2 ,
+            const Array< T , dim > &s3 ,
+            const Array< T , dim > &s4 ,
+            const Array< T , dim > &s5 ,
+            const Array< T , dim > &s6 ,
+            const Array< T , dim > &s7 ,
+            const Array< T , dim > &s8 ,
+            const Array< T , dim > &s9 ,
+            const Array< T , dim > &s10 ,
+            const Array< T , dim > &s11 ,
+            const Array< T , dim > &s12 ,
+            const Array< T , dim > &s13 ,
+            const Array< T , dim > &s14 ,
+            const Array< T , dim > &s15 , Op op )
     {
         for( size_t i=0 ; i<dim ; ++i )
             op( s1[i] , s2[i] , s3[i] , s4[i] , s5[i] , s6[i] , s7[i] , s8[i] , s9[i] , s10[i] , s11[i] , s12[i] , s13[i] , s14[i] , s15[i] );
     }
 
 
-    template< typename T , size_t dim >
-    static typename norm_result_type< boost::array< T , dim > >::type norm_inf( const boost::array< T , dim > &s )
+    template < template < typename, size_t > class Array, typename T,
+               size_t dim>
+    static typename norm_result_type< Array< T , dim > >::type norm_inf( const Array< T , dim > &s )
     {
         BOOST_USING_STD_MAX();
         using std::abs;
-        typedef typename norm_result_type< boost::array< T , dim > >::type result_type;
+        typedef typename norm_result_type< Array< T , dim > >::type result_type;
         result_type init = static_cast< result_type >( 0 );
         for( size_t i=0 ; i<dim ; ++i )
             init = max BOOST_PREVENT_MACRO_SUBSTITUTION ( init , static_cast< result_type >(abs(s[i])) );
diff --git a/boost/numeric/odeint/algebra/multi_array_algebra.hpp b/boost/numeric/odeint/algebra/multi_array_algebra.hpp
index 4d1fc16c3a..0bc476e631 100644
--- a/boost/numeric/odeint/algebra/multi_array_algebra.hpp
+++ b/boost/numeric/odeint/algebra/multi_array_algebra.hpp
@@ -113,13 +113,13 @@ struct multi_array_algebra
     template< class S1 , class S2 , class S3 , class S4 , class S5 , class S6 ,class S7 , class S8 , class S9 , class S10 , class S11 , class S12 , class S13 , class S14 , class Op >
     static void for_each14( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , S5 &s5 , S6 &s6 , S7 &s7 , S8 &s8 , S9 &s9 , S10 &s10 , S11 &s11 , S12 &s12 , S13 &s13 , S14 &s14 , Op op )
     {
-        detail::for_each14( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s14.data() , op );
+        detail::for_each14( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s13.data() , s14.data() , op );
     }
 
     template< class S1 , class S2 , class S3 , class S4 , class S5 , class S6 ,class S7 , class S8 , class S9 , class S10 , class S11 , class S12 , class S13 , class S14 , class S15 , class Op >
     static void for_each15( S1 &s1 , S2 &s2 , S3 &s3 , S4 &s4 , S5 &s5 , S6 &s6 , S7 &s7 , S8 &s8 , S9 &s9 , S10 &s10 , S11 &s11 , S12 &s12 , S13 &s13 , S14 &s14 , S15 &s15 , Op op )
     {
-        detail::for_each15( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s14.data() , s15.data() , op );
+        detail::for_each15( s1.data() , s1.data() + s1.num_elements() , s2.data() , s3.data() , s4.data() , s5.data() , s6.data() , s7.data() , s8.data() , s9.data() , s10.data() , s11.data() , s12.data() , s13.data() , s14.data() , s15.data() , op );
     }
     
     template< typename S >
diff --git a/boost/numeric/odeint/external/eigen/eigen.hpp b/boost/numeric/odeint/external/eigen/eigen.hpp
new file mode 100644
index 0000000000..620acae51f
--- /dev/null
+++ b/boost/numeric/odeint/external/eigen/eigen.hpp
@@ -0,0 +1,27 @@
+/*
+  [auto_generated]
+  boost/numeric/odeint/external/eigen/eigen.hpp
+
+  [begin_description]
+  tba.
+  [end_description]
+
+  Copyright 2009-2012 Karsten Ahnert
+  Copyright 2009-2012 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_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED
+
+
+#include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp>
+#include <boost/numeric/odeint/external/eigen/eigen_resize.hpp>
+
+
+#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_HPP_INCLUDED
diff --git a/boost/numeric/odeint/external/eigen/eigen_algebra.hpp b/boost/numeric/odeint/external/eigen/eigen_algebra.hpp
index 6ca2026b9f..b4ee5c3b1f 100644
--- a/boost/numeric/odeint/external/eigen/eigen_algebra.hpp
+++ b/boost/numeric/odeint/external/eigen/eigen_algebra.hpp
@@ -71,6 +71,19 @@ operator/(const Eigen::MatrixBase<D1> &x1, const Eigen::MatrixBase<D2> &x2) {
     return x1.cwiseQuotient(x2);
 }
 
+
+template< typename D >
+inline const 
+typename Eigen::CwiseUnaryOp<
+    typename Eigen::internal::scalar_abs_op<
+        typename Eigen::internal::traits< D >::Scalar > ,
+        const D >
+abs( const Eigen::MatrixBase< D > &m ) {
+    return m.cwiseAbs();
+}
+
+
+
 } // end Eigen namespace
 
 
diff --git a/boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp b/boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp
new file mode 100644
index 0000000000..6c8a3a276f
--- /dev/null
+++ b/boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp
@@ -0,0 +1,49 @@
+/*
+  [auto_generated]
+  boost/numeric/odeint/external/eigen/eigen_algebra_dispatcher.hpp
+
+  [begin_description]
+  tba.
+  [end_description]
+
+  Copyright 2009-2012 Karsten Ahnert
+  Copyright 2009-2012 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_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED
+
+
+namespace boost {
+namespace numeric {
+namespace odeint {
+
+
+template< class Derived >
+struct algebra_dispatcher_sfinae< Derived ,
+                      typename boost::enable_if< typename boost::is_base_of< Eigen::MatrixBase< Derived > , Derived >::type >::type >
+{ 
+    typedef vector_space_algebra algebra_type;
+};
+
+
+template < class Derived  >
+struct algebra_dispatcher_sfinae< Derived ,
+                      typename boost::enable_if< typename boost::is_base_of< Eigen::ArrayBase< Derived > , Derived >::type >::type >
+{ 
+    typedef vector_space_algebra algebra_type;
+};
+
+
+
+} // namespace odeint
+} // namespace numeric
+} // namespace boost
+
+
+#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_EIGEN_EIGEN_ALGEBRA_DISPATCHER_HPP_INCLUDED
diff --git a/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp b/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp
index 5deba2cb57..0961346350 100644
--- a/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp
+++ b/boost/numeric/odeint/external/thrust/thrust_algebra_dispatcher.hpp
@@ -24,6 +24,7 @@
 #include <boost/numeric/odeint/external/thrust/thrust_algebra.hpp>
 #include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
 
+// specializations for the standard thrust containers
 
 namespace boost {
 namespace numeric {
@@ -48,5 +49,59 @@ struct algebra_dispatcher< thrust::device_vector< T , A > >
 } // namespace boost
 
 
+// add support for thrust backend vectors, if available
+
+#include <thrust/version.h>
+
+#if THRUST_VERSION >= 100600
+
+// specialization for thrust cpp vector
+#include <thrust/system/cpp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct algebra_dispatcher< thrust::cpp::vector< T , A > >
+    {
+        typedef thrust_algebra algebra_type;
+    };
+} } }
+
+// specialization for thrust omp vector
+#ifdef _OPENMP
+#include <thrust/system/omp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct algebra_dispatcher< thrust::omp::vector< T , A > >
+    {
+        typedef thrust_algebra algebra_type;
+    };
+} } }
+#endif // _OPENMP
+
+// specialization for thrust tbb vector
+#ifdef TBB_VERSION_MAJOR
+#include <thrust/system/tbb/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct algebra_dispatcher< thrust::tbb::vector< T , A > >
+    {
+        typedef thrust_algebra algebra_type;
+    };
+} } }
+#endif // TBB_VERSION_MAJOR
+
+// specialization for thrust cuda vector
+#ifdef __CUDACC__
+#include <thrust/system/cuda/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct algebra_dispatcher< thrust::cuda::vector< T , A > >
+    {
+        typedef thrust_algebra algebra_type;
+    };
+} } }
+#endif // __CUDACC__
+
+#endif // THRUST_VERSION >= 100600
+
 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_ALGEBRA_DISPATCHER_HPP_DEFINED
 
diff --git a/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp b/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp
index 0105a4e344..e9b3a64234 100644
--- a/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp
+++ b/boost/numeric/odeint/external/thrust/thrust_operations_dispatcher.hpp
@@ -6,8 +6,8 @@
   operations_dispatcher specialization for thrust
   [end_description]
 
-  Copyright 2013 Karsten Ahnert
-  Copyright 2013 Mario Mulansky
+  Copyright 2013-2014 Karsten Ahnert
+  Copyright 2013-2014 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/external/thrust/thrust_operations.hpp>
 #include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
 
+// support for the standard thrust containers
 
 namespace boost {
 namespace numeric {
@@ -47,6 +48,60 @@ struct operations_dispatcher< thrust::device_vector< T , A > >
 } // namespace numeric
 } // namespace boost
 
+// add support for thrust backend vectors, if available
 
-#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_ALGEBRA_DISPATCHER_HPP_DEFINED
+#include <thrust/version.h>
+
+#if THRUST_VERSION >= 100600
+
+// specialization for thrust cpp vector
+#include <thrust/system/cpp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct operations_dispatcher< thrust::cpp::vector< T , A > >
+    {
+        typedef thrust_operations operations_type;
+    };
+} } }
+
+// specialization for thrust omp vector
+#ifdef _OPENMP
+#include <thrust/system/omp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct operations_dispatcher< thrust::omp::vector< T , A > >
+    {
+        typedef thrust_operations operations_type;
+    };
+} } }
+#endif // _OPENMP
+
+// specialization for thrust tbb vector
+#ifdef TBB_VERSION_MAJOR
+#include <thrust/system/tbb/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct operations_dispatcher< thrust::tbb::vector< T , A > >
+    {
+        typedef thrust_operations operations_type;
+    };
+} } }
+#endif // TBB_VERSION_MAJOR
+
+// specialization for thrust cuda vector
+#ifdef __CUDACC__
+#include <thrust/system/cuda/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    template< class T , class A >
+    struct operations_dispatcher< thrust::cuda::vector< T , A > >
+    {
+        typedef thrust_operations operations_type;
+    };
+} } }
+#endif // __CUDACC__
+
+#endif // THRUST_VERSION >= 100600
+
+
+#endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_OPERATIONS_DISPATCHER_HPP_DEFINED
 
diff --git a/boost/numeric/odeint/external/thrust/thrust_resize.hpp b/boost/numeric/odeint/external/thrust/thrust_resize.hpp
index 8f1b181e25..2f7f7b23b1 100644
--- a/boost/numeric/odeint/external/thrust/thrust_resize.hpp
+++ b/boost/numeric/odeint/external/thrust/thrust_resize.hpp
@@ -6,7 +6,7 @@
  Enable resizing for thrusts device and host_vector.
  [end_description]
 
- Copyright 2010-2012 Mario Mulansky
+ Copyright 2010-2014 Mario Mulansky
  Copyright 2010-2011 Karsten Ahnert
 
  Distributed under the Boost Software License, Version 1.0.
@@ -18,102 +18,168 @@
 #ifndef BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED
 #define BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED
 
+#include <boost/range.hpp>
 
 #include <thrust/device_vector.h>
 #include <thrust/host_vector.h>
+#include <thrust/distance.h>
 
+#include <boost/numeric/odeint/util/resize.hpp>
+#include <boost/numeric/odeint/util/same_size.hpp>
 #include <boost/numeric/odeint/util/copy.hpp>
 
 namespace boost {
 namespace numeric {
 namespace odeint {
 
-template< class T >
-struct is_resizeable< thrust::device_vector< T > >
-{
-    struct type : public boost::true_type { };
-    const static bool value = type::value;
-};
-
-template< class T >
-struct same_size_impl< thrust::device_vector< T > , thrust::device_vector< T > >
-{
-    static bool same_size( const thrust::device_vector< T > &x , const thrust::device_vector< T > &y )
-    {
-        return x.size() == y.size();
-    }
-};
-
-template< class T >
-struct resize_impl< thrust::device_vector< T > , thrust::device_vector< T > >
-{
-    static void resize( thrust::device_vector< T > &x , const thrust::device_vector< T > &y )
-    {
-        x.resize( y.size() );
-    }
-};
-
-
-template< class T >
-struct is_resizeable< thrust::host_vector< T > >
-{
-    struct type : public boost::true_type { };
-    const static bool value = type::value;
-};
-
-template< class T >
-struct same_size_impl< thrust::host_vector< T > , thrust::host_vector< T > >
-{
-    static bool same_size( const thrust::host_vector< T > &x , const thrust::host_vector< T > &y )
-    {
-        return x.size() == y.size();
-    }
-};
-
-template< class T >
-struct resize_impl< thrust::host_vector< T > , thrust::host_vector< T > >
-{
-    static void resize( thrust::host_vector< T > &x , const thrust::host_vector< T > &y )
-    {
-        x.resize( y.size() );
-    }
-};
-
-
-
-template< class Container1, class Value >
-struct copy_impl< Container1 , thrust::device_vector< Value > >
-{
-    static void copy( const Container1 &from , thrust::device_vector< Value > &to )
-    {
-        thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) );
-    }
-};
-
-template< class Value , class Container2 >
-struct copy_impl< thrust::device_vector< Value > , Container2 >
-{
-    static void copy( const thrust::device_vector< Value > &from , Container2 &to )
-    {
-        thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) );
-    }
-};
-
-template< class Value >
-struct copy_impl< thrust::device_vector< Value > , thrust::device_vector< Value > >
-{
-    static void copy( const thrust::device_vector< Value > &from , thrust::device_vector< Value > &to )
-    {
-        thrust::copy( boost::begin( from ) , boost::end( from ) , boost::begin( to ) );
-    }
-};
-
-
+// some macros that define the necessary utilities
+
+#define ODEINT_THRUST_VECTOR_IS_RESIZEABLE( THRUST_VECTOR ) \
+template< class T , class A >                               \
+struct is_resizeable< THRUST_VECTOR<T,A> >                  \
+{                                                           \
+    struct type : public boost::true_type { };              \
+    const static bool value = type::value;                  \
+};                                                          \
+
+#define ODEINT_TRHUST_VECTOR_RESIZE_IMPL( THRUST_VECTOR )     \
+template< class T, class A >                                  \
+struct resize_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> > \
+{                                                             \
+    static void resize( THRUST_VECTOR<T,A> &x ,               \
+                        const THRUST_VECTOR<T,A> &y )         \
+    {                                                         \
+        x.resize( y.size() );                                 \
+    }                                                         \
+};                                                            \
+template< class T, class A, typename Range >                  \
+struct resize_impl< THRUST_VECTOR<T,A> , Range >              \
+{                                                             \
+    static void resize( THRUST_VECTOR<T,A> &x ,               \
+                        const Range &y )                      \
+    {                                                         \
+        x.resize( thrust::distance(boost::begin(y),           \
+                                   boost::end(y)));           \
+    }                                                         \
+};                                                            \
+
+
+#define ODEINT_THRUST_SAME_SIZE_IMPL( THRUST_VECTOR )            \
+template< class T , class A >                                    \
+struct same_size_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> > \
+{                                                                \
+    static bool same_size( const THRUST_VECTOR<T,A> &x ,         \
+                           const THRUST_VECTOR<T,A> &y )         \
+    {                                                            \
+        return x.size() == y.size();                             \
+    }                                                            \
+};                                                               \
+template< class T , class A, typename Range >                    \
+struct same_size_impl< THRUST_VECTOR<T,A> , Range >              \
+{                                                                \
+    static bool same_size( const THRUST_VECTOR<T,A> &x ,         \
+                           const Range &y )                      \
+    {                                                            \
+        return x.size() == thrust::distance(boost::begin(y),     \
+                                            boost::end(y));      \
+    }                                                            \
+};                                                               \
+
+
+#define ODEINT_THRUST_COPY_IMPL( THRUST_VECTOR )                        \
+template< class Container1 , class T , class A >                        \
+struct copy_impl< Container1 , THRUST_VECTOR<T,A> >                     \
+{                                                                       \
+    static void copy( const Container1 &from , THRUST_VECTOR<T,A> &to ) \
+    {                                                                   \
+        thrust::copy( boost::begin( from ) , boost::end( from ) ,       \
+                      boost::begin( to ) );                             \
+    }                                                                   \
+};                                                                      \
+                                                                        \
+template< class T , class A , class Container2 >                        \
+struct copy_impl< THRUST_VECTOR<T,A> , Container2 >                     \
+{                                                                       \
+    static void copy( const THRUST_VECTOR<T,A> &from , Container2 &to ) \
+    {                                                                   \
+        thrust::copy( boost::begin( from ) , boost::end( from ) ,       \
+                      boost::begin( to ) );                             \
+    }                                                                   \
+};                                                                      \
+                                                                        \
+template< class T , class A >                                           \
+struct copy_impl< THRUST_VECTOR<T,A> , THRUST_VECTOR<T,A> >             \
+{                                                                       \
+    static void copy( const THRUST_VECTOR<T,A> &from ,                  \
+                      THRUST_VECTOR<T,A> &to )                          \
+    {                                                                   \
+        thrust::copy( boost::begin( from ) , boost::end( from ) ,       \
+                      boost::begin( to ) );                             \
+    }                                                                   \
+};                                                                      \
+
+// add support for the standard thrust containers
+
+ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::device_vector )
+ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::device_vector )
+ODEINT_THRUST_SAME_SIZE_IMPL( thrust::device_vector )
+ODEINT_THRUST_COPY_IMPL( thrust::device_vector )
+
+
+ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::host_vector )
+ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::host_vector )
+ODEINT_THRUST_SAME_SIZE_IMPL( thrust::host_vector )
+ODEINT_THRUST_COPY_IMPL( thrust::host_vector )
 
 
 } // odeint
 } // numeric
 } // boost
 
+// add support for thrust backend vectors, if available
+
+#include <thrust/version.h>
+
+#if THRUST_VERSION >= 100600
+
+#include <thrust/system/cpp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::cpp::vector )
+    ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::cpp::vector )
+    ODEINT_THRUST_SAME_SIZE_IMPL( thrust::cpp::vector )
+    ODEINT_THRUST_COPY_IMPL( thrust::cpp::vector )
+} } }
+
+#ifdef _OPENMP
+#include <thrust/system/omp/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::omp::vector )
+    ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::omp::vector )
+    ODEINT_THRUST_SAME_SIZE_IMPL( thrust::omp::vector )
+    ODEINT_THRUST_COPY_IMPL( thrust::omp::vector )
+} } }
+#endif // _OPENMP
+
+#ifdef TBB_VERSION_MAJOR
+#include <thrust/system/tbb/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::tbb::vector )
+    ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::tbb::vector )
+    ODEINT_THRUST_SAME_SIZE_IMPL( thrust::tbb::vector )
+    ODEINT_THRUST_COPY_IMPL( thrust::tbb::vector )
+} } }
+#endif // TBB_VERSION_MAJOR
+
+#ifdef __CUDACC__
+#include <thrust/system/cuda/vector.h>
+namespace boost { namespace numeric { namespace odeint {
+    ODEINT_THRUST_VECTOR_IS_RESIZEABLE( thrust::cuda::vector )
+    ODEINT_TRHUST_VECTOR_RESIZE_IMPL( thrust::cuda::vector )
+    ODEINT_THRUST_SAME_SIZE_IMPL( thrust::cuda::vector )
+    ODEINT_THRUST_COPY_IMPL( thrust::cuda::vector )
+} } }
+#endif // __CUDACC__
+
+#endif // THRUST_VERSION >= 100600
 
 #endif // BOOST_NUMERIC_ODEINT_EXTERNAL_THRUST_THRUST_RESIZE_HPP_INCLUDED
diff --git a/boost/numeric/odeint/integrate/integrate.hpp b/boost/numeric/odeint/integrate/integrate.hpp
index 6cb96c0f1e..446656b58a 100644
--- a/boost/numeric/odeint/integrate/integrate.hpp
+++ b/boost/numeric/odeint/integrate/integrate.hpp
@@ -48,6 +48,15 @@ integrate( System system , State &start_state , Time start_time , Time end_time
     return integrate_adaptive( stepper_type() , system , start_state , start_time , end_time , dt , observer );
 }
 
+template< class Value , class System , class State , class Time , class Observer >
+size_t 
+integrate( System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
+{
+    typedef controlled_runge_kutta< runge_kutta_dopri5< State , Value , State , Time > > stepper_type;
+    return integrate_adaptive( stepper_type() , system , start_state , start_time , end_time , dt , observer );
+}
+
+
 
 
 /*
@@ -59,6 +68,13 @@ size_t integrate( System system , State &start_state , Time start_time , Time en
     return integrate( system , start_state , start_time , end_time , dt , null_observer() );
 }
 
+template< class Value , class System , class State , class Time >
+size_t integrate( System system , State &start_state , Time start_time , Time end_time , Time dt )
+{
+    return integrate< Value >( system , start_state , start_time , end_time , dt , null_observer() );
+}
+
+
 
 /**
  * \fn integrate( System system , State &start_state , Time start_time , Time end_time , Time dt , Observer observer )
@@ -70,6 +86,9 @@ size_t integrate( System system , State &start_state , Time start_time , Time en
  * integration with step size control, thus dt changes during the integration.
  * This method uses standard error bounds of 1E-6.
  * After each step, the observer is called.
+ * 
+ * \attention A second version of this function template exists which explicitly
+ * expects the value type as template parameter, i.e. integrate< double >( sys , x , t0 , t1 , dt , obs );
  *
  * \param system The system function to solve, hence the r.h.s. of the 
  * ordinary differential equation.
@@ -92,6 +111,9 @@ size_t integrate( System system , State &start_state , Time start_time , Time en
  * integration with step size control, thus dt changes during the integration.
  * This method uses standard error bounds of 1E-6.
  * No observer is called.
+ * 
+ * \attention A second version of this function template exists which explicitly
+ * expects the value type as template parameter, i.e. integrate< double >( sys , x , t0 , t1 , dt );
  *
  * \param system The system function to solve, hence the r.h.s. of the 
  * ordinary differential equation.
diff --git a/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp b/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp
index 643b53937f..e5b4b3cf48 100644
--- a/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp
+++ b/boost/numeric/odeint/iterator/impl/adaptive_iterator_impl.hpp
@@ -205,7 +205,7 @@ namespace odeint {
          * \param sys The system function (ODE) to solve.
          * \param s The initial state.
          */
-        adaptive_iterator_impl( stepper_type stepper , system_type sys , state_type &s )
+        adaptive_iterator_impl( stepper_type stepper , system_type sys , state_type& /* s */ )
             : base_type( stepper , sys ) { }
 
     protected:
diff --git a/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp b/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp
index 449acfa79f..e23474c7af 100644
--- a/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp
+++ b/boost/numeric/odeint/iterator/impl/const_step_iterator_impl.hpp
@@ -77,7 +77,7 @@ namespace odeint {
          * \param sys The system function (ODE) to solve.
          * \param s The initial state. const_step_iterator stores a reference of s and changes its value during the iteration.
          */
-        const_step_iterator_impl( stepper_type stepper , system_type sys , state_type &s )
+        const_step_iterator_impl( stepper_type stepper , system_type sys , state_type& /* s */ )
             : base_type( stepper , sys ) { }
 
     protected:
diff --git a/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp b/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp
index 147799040d..f90f8875ce 100644
--- a/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp
+++ b/boost/numeric/odeint/iterator/integrate/detail/integrate_n_steps.hpp
@@ -48,7 +48,7 @@ Time integrate_n_steps(
         Observer observer , stepper_tag )
 {
     // ToDo: is there a better way to extract the final time?
-    Time t;
+    Time t = start_time; // Assignment is only here to avoid warnings.
     boost::for_each( make_n_step_time_range( stepper , system , start_state ,
                                              start_time , dt , num_of_steps ) ,
                      obs_caller_time< Observer , Time >( t , observer ) );
@@ -91,7 +91,7 @@ Time integrate_n_steps(
         Observer observer , dense_output_stepper_tag )
 {
     // ToDo: is there a better way to extract the final time?
-    Time t;
+    Time t = start_time;  // Assignment is only here to avoid warnings.
     boost::for_each( make_n_step_time_range( stepper , system , start_state ,
                                              start_time , dt , num_of_steps ) ,
                      obs_caller_time< Observer , Time >( t , observer ) );
diff --git a/boost/numeric/odeint/stepper/adams_bashforth.hpp b/boost/numeric/odeint/stepper/adams_bashforth.hpp
index 59ed1c9025..5ff1e8358c 100644
--- a/boost/numeric/odeint/stepper/adams_bashforth.hpp
+++ b/boost/numeric/odeint/stepper/adams_bashforth.hpp
@@ -37,6 +37,7 @@
 
 #include <boost/numeric/odeint/stepper/stepper_categories.hpp>
 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
+#include <boost/numeric/odeint/stepper/extrapolation_stepper.hpp>
 
 #include <boost/numeric/odeint/stepper/base/algebra_stepper_base.hpp>
 
@@ -44,11 +45,28 @@
 #include <boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp>
 #include <boost/numeric/odeint/stepper/detail/rotating_buffer.hpp>
 
+#include <boost/mpl/arithmetic.hpp>
+#include <boost/mpl/min_max.hpp>
+#include <boost/mpl/equal_to.hpp>
+
+namespace mpl = boost::mpl;
+
 
 namespace boost {
 namespace numeric {
 namespace odeint {
 
+    using mpl::int_;
+
+    /* if N >= 4, returns the smallest even number > N, otherwise returns 4 */
+    template < int N >
+    struct order_helper
+        : mpl::max< typename mpl::eval_if<
+                        mpl::equal_to< mpl::modulus< int_< N >, int_< 2 > >,
+                                       int_< 0 > >,
+                        int_< N >, int_< N + 1 > >::type,
+                    int_< 4 > >::type
+    { };
 
 template<
 size_t Steps ,
@@ -59,7 +77,9 @@ class Time = Value ,
 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
 class Operations = typename operations_dispatcher< State >::operations_type ,
 class Resizer = initially_resizer ,
-class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer >
+class InitializingStepper = extrapolation_stepper< order_helper<Steps>::value, 
+                                                   State, Value, Deriv, Time,
+                                                   Algebra, Operations, Resizer >
 >
 class adams_bashforth : public algebra_stepper_base< Algebra , Operations >
 {
@@ -181,7 +201,8 @@ public :
         {
             if( i != 0 ) m_step_storage.rotate();
             sys( x , m_step_storage[0].m_v , t );
-            stepper.do_step( system , x , m_step_storage[0].m_v , t , dt );
+            stepper.do_step_dxdt_impl( system, x, m_step_storage[0].m_v, t,
+                                       dt );
             t += dt;
         }
         m_steps_initialized = steps;
@@ -222,7 +243,8 @@ private:
         {
             if( m_steps_initialized != 0 ) m_step_storage.rotate();
             sys( in , m_step_storage[0].m_v , t );
-            m_initializing_stepper.do_step( system , in , m_step_storage[0].m_v , t , out , dt );
+            m_initializing_stepper.do_step_dxdt_impl(
+                system, in, m_step_storage[0].m_v, t, out, dt );
             ++m_steps_initialized;
         }
         else
diff --git a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp
index 19b0239990..2f7cc4c6fb 100644
--- a/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp
+++ b/boost/numeric/odeint/stepper/adams_bashforth_moulton.hpp
@@ -49,7 +49,8 @@ class Deriv = State ,
 class Time = Value ,
 class Algebra = typename algebra_dispatcher< State >::algebra_type ,
 class Operations = typename operations_dispatcher< State >::operations_type ,
-class Resizer = initially_resizer
+class Resizer = initially_resizer,
+class InitializingStepper = runge_kutta4< State , Value , Deriv , Time , Algebra , Operations, Resizer >
 >
 class adams_bashforth_moulton
 {
@@ -71,12 +72,13 @@ public :
     typedef Operations operations_type;
     typedef Resizer resizer_type;
     typedef stepper_tag stepper_category;
+    typedef InitializingStepper initializing_stepper_type;
 
     static const size_t steps = Steps;
 #ifndef DOXYGEN_SKIP
-    typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_bashforth_type;
+    typedef adams_bashforth< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type, initializing_stepper_type > adams_bashforth_type;
     typedef adams_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > adams_moulton_type;
-    typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > stepper_type;
+    typedef adams_bashforth_moulton< steps , state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type , initializing_stepper_type> stepper_type;
 #endif //DOXYGEN_SKIP
     typedef unsigned short order_type;
     static const order_type order_value = steps;
@@ -158,7 +160,7 @@ private:
         {
             m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl< StateInOut > , detail::ref( *this ) , detail::_1 ) );
             m_adams_bashforth.do_step( system , x , t , m_x.m_v , dt );
-            m_adams_moulton.do_step( system , x , m_x.m_v , t , x , dt , m_adams_bashforth.step_storage() );
+            m_adams_moulton.do_step( system , x , m_x.m_v , t+dt , x , dt , m_adams_bashforth.step_storage() );
         }
         else
         {
diff --git a/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp b/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp
index 3ba11a4fe0..08009dc16b 100644
--- a/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp
+++ b/boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp
@@ -152,6 +152,22 @@ public:
 
 
     /*
+     * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version that is not disabled. Meant for internal use only.
+     */
+    template < class System, class StateInOut, class DerivIn >
+    void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
+                            time_type t, time_type dt )
+    {
+        this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
+    }
+
+
+
+    /*
      * Version 3 : do_step( sys , in , t , out , dt )
      *
      * this version does not solve the forwarding problem, boost.range can not be used
@@ -181,10 +197,21 @@ public:
     {
         this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
     }
-
-
-
-
+    
+    /*
+     * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version that is not disabled. Meant for internal use only.
+     */
+    template < class System, class StateIn, class DerivIn, class StateOut >
+    void do_step_dxdt_impl( System system, const StateIn &in,
+                            const DerivIn &dxdt, time_type t, StateOut &out,
+                            time_type dt )
+    {
+        this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
+    }
 
     /*
      * Version  5 :do_step( sys , x , t , dt , xerr )
diff --git a/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp b/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp
index a055c7fa5f..b1d751a0ce 100644
--- a/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp
+++ b/boost/numeric/odeint/stepper/base/explicit_error_stepper_fsal_base.hpp
@@ -151,11 +151,27 @@ public:
 
 
     /*
+     * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version that is not disabled. Meant for internal use only.
+     */
+    template< class System , class StateInOut , class DerivInOut >
+    void do_step_dxdt_impl( System system , StateInOut &x , DerivInOut &dxdt , time_type t , time_type dt )
+    {
+        m_first_call = true;
+        this->stepper().do_step_impl( system , x , dxdt , t , x , dxdt , dt );
+    }
+
+    /*
      * version 3 : do_step( sys , in , t , out , dt )
      *
-     * this version does not solve the forwarding problem, boost.range can not be used
+     * this version does not solve the forwarding problem, boost.range can not
+     * be used.
      *
-     * the disable is needed to avoid ambiguous overloads if state_type = time_type
+     * the disable is needed to avoid ambiguous overloads if 
+     * state_type = time_type
      */
     template< class System , class StateIn , class StateOut >
     typename boost::disable_if< boost::is_same< StateIn , time_type > , void >::type
@@ -174,12 +190,14 @@ public:
      *
      * this version does not solve the forwarding problem, boost.range can not be used
      */
-    template< class System , class StateIn , class DerivIn , class StateOut , class DerivOut >
-    void do_step( System system , const StateIn &in , const DerivIn &dxdt_in , time_type t ,
-            StateOut &out , DerivOut &dxdt_out , time_type dt )
+    template< class System, class StateIn, class DerivIn, class StateOut,
+               class DerivOut >
+    void do_step( System system, const StateIn &in, const DerivIn &dxdt_in,
+                  time_type t, StateOut &out, DerivOut &dxdt_out, time_type dt )
     {
         m_first_call = true;
-        this->stepper().do_step_impl( system , in , dxdt_in , t , out , dxdt_out , dt );
+        this->stepper().do_step_impl( system, in, dxdt_in, t, out, dxdt_out,
+                                      dt );
     }
 
 
diff --git a/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp b/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp
index 40aab80399..d81c8c7a71 100644
--- a/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp
+++ b/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp
@@ -139,6 +139,21 @@ public:
 
 
     /*
+     * named Version 2: do_step_dxdt_impl( sys , in , dxdt , t , dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version that is not disabled. Meant for internal use only.
+     */
+    template < class System, class StateInOut, class DerivIn >
+    void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
+                            time_type t, time_type dt )
+    {
+        this->stepper().do_step_impl( system , x , dxdt , t , x , dt );
+    }
+
+
+    /*
      * Version 3 : do_step( sys , in , t , out , dt )
      *
      * this version does not solve the forwarding problem, boost.range can not be used
@@ -164,6 +179,22 @@ public:
         this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
     }
 
+
+    /*
+     * named Version 4: do_step_dxdt_impl( sys , in , dxdt , t , out, dt )
+     *
+     * this version is needed when this stepper is used for initializing 
+     * multistep stepper like adams-bashforth. Hence we provide an explicitely
+     * named version. Meant for internal use only.
+     */
+    template < class System, class StateIn, class DerivIn, class StateOut >
+    void do_step_dxdt_impl( System system, const StateIn &in,
+                            const DerivIn &dxdt, time_type t, StateOut &out,
+                            time_type dt )
+    {
+        this->stepper().do_step_impl( system , in , dxdt , t , out , dt );
+    }
+
     template< class StateIn >
     void adjust_size( const StateIn &x )
     {
diff --git a/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp b/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp
index ea3523e5bb..eb09aefc8f 100644
--- a/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp
+++ b/boost/numeric/odeint/stepper/base/symplectic_rkn_stepper_base.hpp
@@ -182,7 +182,7 @@ private:
 
     // stepper for systems with function for dq/dt = f(p) and dp/dt = -f(q)
     template< class System , class StateIn , class StateOut >
-    void do_step_impl( System system , const StateIn &in , time_type t , StateOut &out , time_type dt , boost::mpl::true_ )
+    void do_step_impl( System system , const StateIn &in , time_type /* t */ , StateOut &out , time_type dt , boost::mpl::true_ )
     {
         typedef typename odeint::unwrap_reference< System >::type system_type;
         typedef typename odeint::unwrap_reference< typename system_type::first_type >::type coor_deriv_func_type;
diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
index b604f86a27..e71008a993 100644
--- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
+++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp
@@ -194,7 +194,6 @@ public:
 
         static const value_type val1( 1.0 );
 
-        typename odeint::unwrap_reference< System >::type &sys = system;
         if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) )
         {
             reset(); // system resized -> reset
@@ -219,12 +218,12 @@ public:
             m_midpoint.set_steps( m_interval_sequence[k] );
             if( k == 0 )
             {
-                m_midpoint.do_step( sys , in , dxdt , t , out , dt );
+                m_midpoint.do_step( system , in , dxdt , t , out , dt );
                 /* the first step, nothing more to do */
             }
             else
             {
-                m_midpoint.do_step( sys , in , dxdt , t , m_table[k-1].m_v , dt );
+                m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt );
                 extrapolate( k , m_table , m_coeff , out );
                 // get error estimate
                 m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v ,
@@ -341,7 +340,7 @@ public:
         resize_m_dxdt( x );
         resize_m_xnew( x );
         resize_impl( x );
-        m_midpoint.adjust_size();
+        m_midpoint.adjust_size( x );
     }
 
 
diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
index c352c9605a..0036d1d2a6 100644
--- a/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
+++ b/boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp
@@ -155,8 +155,6 @@ public:
         
         static const value_type val1( 1.0 );
 
-        typename odeint::unwrap_reference< System >::type &sys = system;
-
         bool reject( true );
 
         time_vector h_opt( m_k_max+1 );
@@ -261,6 +259,7 @@ public:
         {
 
             //calculate dxdt for next step and dense output
+            typename odeint::unwrap_reference< System >::type &sys = system;
             sys( out , dxdt_new , t+dt );
 
             //prepare dense output
@@ -367,7 +366,7 @@ public:
     void adjust_size( const StateIn &x )
     {
         resize_impl( x );
-        m_midpoint.adjust_size();
+        m_midpoint.adjust_size( x );
     }
 
 
diff --git a/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp b/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp
index 168bc69202..cc279abbe0 100644
--- a/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp
+++ b/boost/numeric/odeint/stepper/detail/adams_bashforth_call_algebra.hpp
@@ -117,7 +117,7 @@ struct adams_bashforth_call_algebra< 7 , Algebra , Operations >
     {
         //BOOST_ASSERT( false ); // not implemented
         typedef typename Coefficients::value_type value_type;
-        Algebra::for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v ,
+        algebra.for_each9( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v ,
                             typename Operations::template scale_sum8< value_type , Time , Time , Time , Time , Time , Time >(
                                 1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] ) );
     }
@@ -132,7 +132,7 @@ struct adams_bashforth_call_algebra< 8 , Algebra , Operations >
     {
         //BOOST_ASSERT( false ); // not implemented
         typedef typename Coefficients::value_type value_type;
-        Algebra::for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v ,
+        algebra.for_each10( out , in , steps[0].m_v , steps[1].m_v , steps[2].m_v , steps[3].m_v , steps[4].m_v , steps[5].m_v , steps[6].m_v , steps[7].m_v ,
                              typename Operations::template scale_sum9< value_type , Time , Time , Time , Time , Time , Time , Time >(
                                  1.0 , dt * coef[0] , dt * coef[1] , dt * coef[2] , dt * coef[3] , dt * coef[4] , dt * coef[5] , dt * coef[6] , dt * coef[7] ) );
     }
diff --git a/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp b/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp
index f8a1505510..b6f5f2a40b 100644
--- a/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp
+++ b/boost/numeric/odeint/stepper/detail/adams_moulton_call_algebra.hpp
@@ -32,7 +32,7 @@ template< class Algebra , class Operations >
 struct adams_moulton_call_algebra< 1 , Algebra , Operations >
 {
     template< class StateIn , class StateOut , class DerivIn , class StepStorage , class Coefficients , class Time >
-    void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage &steps , const Coefficients &coef , Time dt ) const
+    void operator()( Algebra &algebra , const StateIn &in , StateOut &out , const DerivIn &dxdt , const StepStorage& /* steps */ , const Coefficients &coef , Time dt ) const
     {
         typedef typename Coefficients::value_type value_type;
         algebra.for_each3( out , in , dxdt , typename Operations::template scale_sum2< value_type , Time >( 1.0 , dt * coef[0] ) );
diff --git a/boost/numeric/odeint/stepper/euler.hpp b/boost/numeric/odeint/stepper/euler.hpp
index 443f942763..1c7c126b96 100644
--- a/boost/numeric/odeint/stepper/euler.hpp
+++ b/boost/numeric/odeint/stepper/euler.hpp
@@ -76,7 +76,7 @@ public :
     { }
 
     template< class System , class StateIn , class DerivIn , class StateOut >
-    void do_step_impl( System system , const StateIn &in , const DerivIn &dxdt , time_type t , StateOut &out , time_type dt )
+    void do_step_impl( System /* system */ , const StateIn &in , const DerivIn &dxdt , time_type /* t */ , StateOut &out , time_type dt )
     {
         stepper_base_type::m_algebra.for_each3( out , in , dxdt ,
                 typename operations_type::template scale_sum2< value_type , time_type >( 1.0 , dt ) );
@@ -84,7 +84,7 @@ public :
     }
 
     template< class StateOut , class StateIn1 , class StateIn2 >
-    void calc_state( StateOut &x , time_type t ,  const StateIn1 &old_state , time_type t_old , const StateIn2 &current_state , time_type t_new ) const
+    void calc_state( StateOut &x , time_type t ,  const StateIn1 &old_state , time_type t_old , const StateIn2 & /*current_state*/ , time_type /* t_new */ ) const
     {
         const time_type delta = t - t_old;
         stepper_base_type::m_algebra.for_each3( x , old_state , stepper_base_type::m_dxdt.m_v ,
diff --git a/boost/numeric/odeint/stepper/extrapolation_stepper.hpp b/boost/numeric/odeint/stepper/extrapolation_stepper.hpp
new file mode 100644
index 0000000000..3b688e958b
--- /dev/null
+++ b/boost/numeric/odeint/stepper/extrapolation_stepper.hpp
@@ -0,0 +1,293 @@
+/*
+  [auto_generated]
+  boost/numeric/odeint/stepper/extrapolation_stepper.hpp
+
+  [begin_description]
+  extrapolation stepper
+  [end_description]
+
+  Copyright 2009-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_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED
+#define BOOST_NUMERIC_ODEINT_STEPPER_EXTRAPOLATION_STEPPER_HPP_INCLUDED
+
+#include <iostream>
+
+#include <algorithm>
+
+#include <boost/config.hpp> // for min/max guidelines
+#include <boost/static_assert.hpp>
+
+#include <boost/numeric/odeint/util/bind.hpp>
+#include <boost/numeric/odeint/util/unwrap_reference.hpp>
+
+#include <boost/numeric/odeint/stepper/base/explicit_error_stepper_base.hpp>
+#include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
+#include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
+#include <boost/numeric/odeint/algebra/range_algebra.hpp>
+#include <boost/numeric/odeint/algebra/default_operations.hpp>
+#include <boost/numeric/odeint/algebra/algebra_dispatcher.hpp>
+#include <boost/numeric/odeint/algebra/operations_dispatcher.hpp>
+
+#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/unit_helper.hpp>
+#include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
+
+namespace boost
+{
+namespace numeric
+{
+namespace odeint
+{
+
+template < unsigned short Order, class State, class Value = double,
+           class Deriv = State, class Time = Value,
+           class Algebra = typename algebra_dispatcher< State >::algebra_type,
+           class Operations =
+               typename operations_dispatcher< State >::operations_type,
+           class Resizer = initially_resizer >
+#ifndef DOXYGEN_SKIP
+class extrapolation_stepper
+    : public explicit_error_stepper_base<
+          extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra,
+                                 Operations, Resizer >,
+          Order, Order, Order - 2, State, Value, Deriv, Time, Algebra,
+          Operations, Resizer >
+#else
+class extrapolation_stepper : public explicit_error_stepper_base
+#endif
+{
+
+  private:
+    // check for Order being odd
+    BOOST_STATIC_ASSERT_MSG(
+        ( ( Order % 2 ) == 0 ) && ( Order > 2 ),
+        "extrapolation_stepper requires even Order larger than 2" );
+
+  public:
+#ifndef DOXYGEN_SKIP
+    typedef explicit_error_stepper_base<
+        extrapolation_stepper< Order, State, Value, Deriv, Time, Algebra,
+                               Operations, Resizer >,
+        Order, Order, Order - 2, State, Value, Deriv, Time, Algebra, Operations,
+        Resizer > stepper_base_type;
+#else
+    typedef explicit_error_stepper_base< extrapolation_stepper< ... >, ... >
+    stepper_base_type;
+#endif
+
+    typedef typename stepper_base_type::state_type state_type;
+    typedef typename stepper_base_type::value_type value_type;
+    typedef typename stepper_base_type::deriv_type deriv_type;
+    typedef typename stepper_base_type::time_type time_type;
+    typedef typename stepper_base_type::algebra_type algebra_type;
+    typedef typename stepper_base_type::operations_type operations_type;
+    typedef typename stepper_base_type::resizer_type resizer_type;
+
+#ifndef DOXYGEN_SKIP
+    typedef typename stepper_base_type::stepper_type stepper_type;
+    typedef typename stepper_base_type::wrapped_state_type wrapped_state_type;
+    typedef typename stepper_base_type::wrapped_deriv_type wrapped_deriv_type;
+
+    typedef std::vector< value_type > value_vector;
+    typedef std::vector< value_vector > value_matrix;
+    typedef std::vector< size_t > int_vector;
+    typedef std::vector< wrapped_state_type > state_table_type;
+    typedef modified_midpoint< state_type, value_type, deriv_type, time_type,
+                               algebra_type, operations_type,
+                               resizer_type > midpoint_stepper_type;
+
+#endif // DOXYGEN_SKIP
+
+    typedef unsigned short order_type;
+    static const order_type order_value = stepper_base_type::order_value;
+    static const order_type stepper_order_value =
+        stepper_base_type::stepper_order_value;
+    static const order_type error_order_value =
+        stepper_base_type::error_order_value;
+
+    const static size_t m_k_max = ( order_value - 2 ) / 2;
+
+    extrapolation_stepper( const algebra_type &algebra = algebra_type() )
+        : stepper_base_type( algebra ), m_interval_sequence( m_k_max + 1 ),
+          m_coeff( m_k_max + 1 ), m_table( m_k_max )
+    {
+        for ( unsigned short i = 0; i < m_k_max + 1; i++ )
+        {
+            m_interval_sequence[i] = 2 * ( i + 1 );
+            m_coeff[i].resize( i );
+            for ( size_t k = 0; k < i; ++k )
+            {
+                const value_type r =
+                    static_cast< value_type >( m_interval_sequence[i] ) /
+                    static_cast< value_type >( m_interval_sequence[k] );
+                m_coeff[i][k] =
+                    static_cast< value_type >( 1 ) /
+                    ( r * r - static_cast< value_type >(
+                                  1 ) ); // coefficients for extrapolation
+            }
+        }
+    }
+
+    template < class System, class StateIn, class DerivIn, class StateOut,
+               class Err >
+    void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt,
+                       time_type t, StateOut &out, time_type dt, Err &xerr )
+    {
+        // std::cout << "dt: " << dt << std::endl;
+        // normal step
+        do_step_impl( system, in, dxdt, t, out, dt );
+
+        static const value_type val1( 1.0 );
+        // additionally, perform the error calculation
+        stepper_base_type::m_algebra.for_each3(
+            xerr, out, m_table[0].m_v,
+            typename operations_type::template scale_sum2<
+                value_type, value_type >( val1, -val1 ) );
+    }
+
+    template < class System, class StateInOut, class DerivIn, class Err >
+    void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
+                          time_type t, time_type dt, Err &xerr )
+    {
+        // normal step
+        do_step_impl_io( system, inout, dxdt, t, dt );
+
+        static const value_type val1( 1.0 );
+        // additionally, perform the error calculation
+        stepper_base_type::m_algebra.for_each3(
+            xerr, inout, m_table[0].m_v,
+            typename operations_type::template scale_sum2<
+                value_type, value_type >( val1, -val1 ) );
+    }
+
+    template < class System, class StateIn, class DerivIn, class StateOut >
+    void do_step_impl( System system, const StateIn &in, const DerivIn &dxdt,
+                       time_type t, StateOut &out, time_type dt )
+    {
+        m_resizer.adjust_size(
+            in, detail::bind( &stepper_type::template resize_impl< StateIn >,
+                              detail::ref( *this ), detail::_1 ) );
+        size_t k = 0;
+        m_midpoint.set_steps( m_interval_sequence[k] );
+        m_midpoint.do_step( system, in, dxdt, t, out, dt );
+        for ( k = 1; k <= m_k_max; ++k )
+        {
+            m_midpoint.set_steps( m_interval_sequence[k] );
+            m_midpoint.do_step( system, in, dxdt, t, m_table[k - 1].m_v, dt );
+            extrapolate( k, m_table, m_coeff, out );
+        }
+    }
+
+    template < class System, class StateInOut, class DerivIn >
+    void do_step_impl_io( System system, StateInOut &inout, const DerivIn &dxdt,
+                          time_type t, time_type dt )
+    {
+        // special care for inout
+        m_xout_resizer.adjust_size(
+            inout,
+            detail::bind( &stepper_type::template resize_m_xout< StateInOut >,
+                          detail::ref( *this ), detail::_1 ) );
+        do_step_impl( system, inout, dxdt, t, m_xout.m_v, dt );
+        boost::numeric::odeint::copy( m_xout.m_v, inout );
+    }
+
+    template < class System, class StateInOut, class DerivIn >
+    void do_step_dxdt_impl( System system, StateInOut &x, const DerivIn &dxdt,
+                            time_type t, time_type dt )
+    {
+        do_step_impl_io( system , x , dxdt , t , dt );
+    }
+
+    template < class System, class StateIn, class DerivIn, class StateOut >
+    void do_step_dxdt_impl( System system, const StateIn &in,
+                            const DerivIn &dxdt, time_type t, StateOut &out,
+                            time_type dt )
+    {
+        do_step_impl( system , in , dxdt , t , out , dt );
+    }
+
+
+    template < class StateIn > void adjust_size( const StateIn &x )
+    {
+        resize_impl( x );
+        m_midpoint.adjust_size( x );
+    }
+
+  private:
+    template < class StateIn > bool resize_impl( const StateIn &x )
+    {
+        bool resized( false );
+        for ( size_t i = 0; i < m_k_max; ++i )
+            resized |= adjust_size_by_resizeability(
+                m_table[i], x, typename is_resizeable< state_type >::type() );
+        return resized;
+    }
+
+    template < class StateIn > bool resize_m_xout( const StateIn &x )
+    {
+        return adjust_size_by_resizeability(
+            m_xout, x, typename is_resizeable< state_type >::type() );
+    }
+
+    template < class StateInOut >
+    void extrapolate( size_t k, state_table_type &table,
+                      const value_matrix &coeff, StateInOut &xest )
+    /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf
+       uses the obtained intermediate results to extrapolate to dt->0
+    */
+    {
+        static const value_type val1 = static_cast< value_type >( 1.0 );
+
+        for ( int j = k - 1; j > 0; --j )
+        {
+            stepper_base_type::m_algebra.for_each3(
+                table[j - 1].m_v, table[j].m_v, table[j - 1].m_v,
+                typename operations_type::template scale_sum2<
+                    value_type, value_type >( val1 + coeff[k][j],
+                                              -coeff[k][j] ) );
+        }
+        stepper_base_type::m_algebra.for_each3(
+            xest, table[0].m_v, xest,
+            typename operations_type::template scale_sum2<
+                value_type, value_type >( val1 + coeff[k][0], -coeff[k][0] ) );
+    }
+
+  private:
+    midpoint_stepper_type m_midpoint;
+
+    resizer_type m_resizer;
+    resizer_type m_xout_resizer;
+
+    int_vector m_interval_sequence; // stores the successive interval counts
+    value_matrix m_coeff;
+
+    wrapped_state_type m_xout;
+    state_table_type m_table; // sequence of states for extrapolation
+};
+
+/******** DOXYGEN *******/
+
+/**
+ * \class extrapolation_stepper
+ * \brief Extrapolation stepper with configurable order, and error estimation.
+ *
+ * The extrapolation stepper is a stepper with error estimation and configurable
+ * order. The order is given as template parameter and needs to be an _odd_
+ * number. The stepper is based on several executions of the modified midpoint
+ * method and a Richardson extrapolation. This is essentially the same technique
+ * as for bulirsch_stoer, but without the variable order.
+ *
+ * \note The Order parameter has to be an even number greater 2.
+ */
+}
+}
+}
+#endif
diff --git a/boost/numeric/odeint/stepper/implicit_euler.hpp b/boost/numeric/odeint/stepper/implicit_euler.hpp
index f3019889cb..e1c64164bc 100644
--- a/boost/numeric/odeint/stepper/implicit_euler.hpp
+++ b/boost/numeric/odeint/stepper/implicit_euler.hpp
@@ -144,7 +144,7 @@ private:
     void solve( state_type &x , matrix_type &m )
     {
         int res = boost::numeric::ublas::lu_factorize( m , m_pm.m_v );
-        if( res != 0 ) exit(0);
+        if( res != 0 ) std::exit(0);
         boost::numeric::ublas::lu_substitute( m , m_pm.m_v , x );
     }
 
diff --git a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
index 4b3824fc08..df4e6c48b6 100644
--- a/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
+++ b/boost/numeric/odeint/stepper/rosenbrock4_controller.hpp
@@ -66,6 +66,7 @@ public:
     {
         BOOST_USING_STD_MAX();
         using std::abs;
+        using std::sqrt;
         
         const size_t n = x.size();
         value_type err = 0.0 , sk = 0.0;
diff --git a/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp b/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp
index 2b8e70b166..32bda0bd7f 100644
--- a/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp
+++ b/boost/numeric/odeint/stepper/runge_kutta4_classic.hpp
@@ -121,11 +121,20 @@ public :
 
         // dt * m_dxh = k4
         sys( m_x_tmp.m_v , m_dxh.m_v , t + dt );
+
         //x += dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
         time_type dt6 = dt / static_cast< value_type >( 6 );
         time_type dt3 = dt / static_cast< value_type >( 3 );
         stepper_base_type::m_algebra.for_each6( out , in , dxdt , m_dxt.m_v , m_dxm.m_v , m_dxh.m_v ,
-                typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) );
+                                             typename operations_type::template scale_sum5< value_type , time_type , time_type , time_type , time_type >( 1.0 , dt6 , dt3 , dt3 , dt6 ) );
+        
+        // x += dt/6 * m_dxdt + dt/3 * m_dxt )
+        // stepper_base_type::m_algebra.for_each4( out , in , dxdt , m_dxt.m_v , 
+        //                                         typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt6 , dt3 ) ); 
+        // // x += dt/3 * m_dxm + dt/6 * m_dxh )
+        // stepper_base_type::m_algebra.for_each4( out , out , m_dxm.m_v , m_dxh.m_v , 
+        //                                         typename operations_type::template scale_sum3< value_type , time_type , time_type >( 1.0 , dt3 , dt6 ) ); 
+
     }
 
     template< class StateType >
diff --git a/boost/numeric/odeint/stepper/velocity_verlet.hpp b/boost/numeric/odeint/stepper/velocity_verlet.hpp
index 24f3c0e3b6..3a20fc25c9 100644
--- a/boost/numeric/odeint/stepper/velocity_verlet.hpp
+++ b/boost/numeric/odeint/stepper/velocity_verlet.hpp
@@ -114,7 +114,7 @@ public:
 
         algebra_stepper_base_type::m_algebra.for_each4(
             qout , qin , pin , ain ,
-            typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one_half * dt , one * dt * dt ) );
+            typename operations_type::template scale_sum3< value_type , time_type , time_square_type >( one , one * dt , one_half * dt * dt ) );
 
         typename odeint::unwrap_reference< System >::type & sys = system;
 
diff --git a/boost/numeric/odeint/util/bind.hpp b/boost/numeric/odeint/util/bind.hpp
index 966fd0397b..1201afab84 100644
--- a/boost/numeric/odeint/util/bind.hpp
+++ b/boost/numeric/odeint/util/bind.hpp
@@ -41,9 +41,12 @@ using namespace ::std::placeholders;
 
 #else
 
+// unnamed namespace to avoid multiple declarations (#138)
+namespace {
 using ::boost::bind;
 boost::arg<1> _1;
 boost::arg<2> _2;
+}
 // using ::boost::bind;
 // using ::_1;
 // using ::_2;
diff --git a/boost/numeric/odeint/util/resize.hpp b/boost/numeric/odeint/util/resize.hpp
index db2ce208b0..75e9eff448 100644
--- a/boost/numeric/odeint/util/resize.hpp
+++ b/boost/numeric/odeint/util/resize.hpp
@@ -33,8 +33,8 @@
 namespace boost {
 namespace numeric {
 namespace odeint {
-    
-    
+
+
 template< class StateOut , class StateIn , class Enabler = void >
 struct resize_impl_sfinae
 {
diff --git a/boost/numeric/odeint/util/same_instance.hpp b/boost/numeric/odeint/util/same_instance.hpp
index dc55db30a3..a889ee1925 100644
--- a/boost/numeric/odeint/util/same_instance.hpp
+++ b/boost/numeric/odeint/util/same_instance.hpp
@@ -25,7 +25,7 @@ namespace odeint {
 template< class T1 , class T2 , class Enabler=void >
 struct same_instance_impl
 { 
-    static bool same_instance( const T1 &x1 , const T2 &x2 )
+    static bool same_instance( const T1& /* x1 */ , const T2& /* x2 */ )
     {
         return false;
     }
diff --git a/boost/numeric/odeint/util/unwrap_reference.hpp b/boost/numeric/odeint/util/unwrap_reference.hpp
index d622a5c007..bc7d423a09 100644
--- a/boost/numeric/odeint/util/unwrap_reference.hpp
+++ b/boost/numeric/odeint/util/unwrap_reference.hpp
@@ -31,9 +31,9 @@
 namespace boost {
 
 #if BOOST_NUMERIC_ODEINT_CXX11
-template<typename T> class reference_wrapper;
+template<typename T> struct reference_wrapper;
 
-template<typename T> class unwrap_reference;
+template<typename T> struct unwrap_reference;
 #endif
 
 namespace numeric {
@@ -43,24 +43,21 @@ namespace odeint {
 #if BOOST_NUMERIC_ODEINT_CXX11
 
 template<typename T>
-class unwrap_reference
+struct unwrap_reference
 {
-public:
     typedef typename std::remove_reference<T>::type type;
 };
 
 template<typename T>
-class unwrap_reference< std::reference_wrapper<T> >
+struct unwrap_reference< std::reference_wrapper<T> >
 {
-public:
     typedef typename std::remove_reference<T>::type type;
 };
 
 template<typename T>
-class unwrap_reference< boost::reference_wrapper<T> >
+struct unwrap_reference< boost::reference_wrapper<T> >
 {
-public:
-        typedef typename boost::unwrap_reference<T>::type type;
+    typedef typename boost::unwrap_reference<T>::type type;
 };
 
 #else