diff --git a/examples/multiprecision/Jamfile b/examples/multiprecision/Jamfile
new file mode 100644
index 0000000..9708d43
--- /dev/null
+++ b/examples/multiprecision/Jamfile
@@ -0,0 +1,16 @@
+# Copyright 2011-2013 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# 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)
+
+
+project
+    : requirements 
+      <define>BOOST_ALL_NO_LIB=1
+    : 
+    ;
+    
+
+exe lorenz_mp : lorenz_mp.cpp ;
+exe cmp_precision : cmp_precision.cpp ;
diff --git a/examples/multiprecision/cmp_precision.cpp b/examples/multiprecision/cmp_precision.cpp
new file mode 100644
index 0000000..779b5e2
--- /dev/null
+++ b/examples/multiprecision/cmp_precision.cpp
@@ -0,0 +1,68 @@
+/* Boost libs/numeric/odeint/examples/multiprecision/cmp_precision.cpp
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+
+ example comparing double to multiprecision using Boost.Multiprecision
+
+ 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)
+ */
+
+
+#include <iostream>
+#include <boost/numeric/odeint.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef boost::multiprecision::cpp_dec_float_50 mp_50;
+
+/* we solve the simple ODE x' = 3/(2t^2) + x/(2t)
+ * with initial condition x(1) = 0.
+ * Analytic solution is x(t) = sqrt(t) - 1/t
+ */
+
+void rhs_m( const mp_50 x , mp_50 &dxdt , const mp_50 t )
+{   // version for multiprecision
+    dxdt = mp_50(3)/(mp_50(2)*t*t) + x/(mp_50(2)*t);
+}
+
+void rhs_d( const double x , double &dxdt , const double t )
+{   // version for double precision
+    dxdt = 3.0/(2.0*t*t) + x/(2.0*t);
+}
+
+// state_type = mp_50 = deriv_type = time_type = mp_50
+typedef runge_kutta4< mp_50 , mp_50 , mp_50 , mp_50 , vector_space_algebra , default_operations , never_resizer > stepper_type_m;
+
+typedef runge_kutta4< double , double , double , double , vector_space_algebra , default_operations , never_resizer > stepper_type_d;
+
+int main()
+{
+
+    stepper_type_m stepper_m;
+    stepper_type_d stepper_d;
+
+    mp_50 dt_m( 0.5 );
+    double dt_d( 0.5 );
+
+    cout << "dt" << '\t' << "mp" << '\t' << "double" << endl;
+    
+    while( dt_m > 1E-20 )
+    {
+
+        mp_50 x_m = 0; //initial value x(1) = 0
+        stepper_m.do_step( rhs_m , x_m , mp_50( 1 ) , dt_m );
+        double x_d = 0;
+        stepper_d.do_step( rhs_d , x_d , 1.0 , dt_d );        
+
+        cout << dt_m << '\t';
+        cout << abs((x_m - (sqrt(1+dt_m)-mp_50(1)/(1+dt_m)))/x_m) << '\t' ;
+        cout << abs((x_d - (sqrt(1+dt_d)-mp_50(1)/(1+dt_d)))/x_d) << endl ;
+        dt_m /= 2;
+        dt_d /= 2;
+    }
+}
diff --git a/examples/multiprecision/lorenz_mp.cpp b/examples/multiprecision/lorenz_mp.cpp
new file mode 100644
index 0000000..c3e5b82
--- /dev/null
+++ b/examples/multiprecision/lorenz_mp.cpp
@@ -0,0 +1,81 @@
+/*
+ * lorenz_mp.cpp
+ *
+ * This example demonstrates how odeint can be used with boost.multiprecision.
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2013 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)
+ */
+
+
+
+#include <iostream>
+
+//[ mp_lorenz_defs
+#include <boost/numeric/odeint.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef boost::multiprecision::cpp_dec_float_50 value_type;
+
+typedef boost::array< value_type , 3 > state_type;
+//]
+
+//[ mp_lorenz_rhs
+struct lorenz
+{
+    void operator()( const state_type &x , state_type &dxdt , value_type t ) const
+    {
+        const value_type sigma( 10 );
+        const value_type R( 28 );
+        const value_type b( value_type( 8 ) / value_type( 3 ) );
+
+        dxdt[0] = sigma * ( x[1] - x[0] );
+        dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+        dxdt[2] = -b * x[2] + x[0] * x[1];
+    }
+};
+//]
+
+
+
+
+struct streaming_observer
+{
+    std::ostream& m_out;
+
+    streaming_observer( std::ostream &out ) : m_out( out ) { }
+
+    template< class State , class Time >
+    void operator()( const State &x , Time t ) const
+    {
+        m_out << t;
+        for( size_t i=0 ; i<x.size() ; ++i ) m_out << "\t" << x[i] ;
+        m_out << "\n";
+    }
+};
+
+
+
+
+
+
+int main( int argc , char **argv )
+{
+    //[ mp_lorenz_int
+    state_type x = {{ value_type( 10.0 ) , value_type( 10.0 ) , value_type( 10.0 ) }};
+
+    cout.precision( 50 );
+    integrate_const( runge_kutta4< state_type , value_type >() ,
+            lorenz() , x , value_type( 0.0 ) , value_type( 10.0 ) , value_type( value_type( 1.0 ) / value_type( 10.0 ) ) ,
+            streaming_observer( cout ) );
+    //]
+
+    return 0;
+}