diff --git a/examples/openmp/Jamfile.v2 b/examples/openmp/Jamfile.v2
new file mode 100644
index 0000000..cef4f67
--- /dev/null
+++ b/examples/openmp/Jamfile.v2
@@ -0,0 +1,23 @@
+# Copyright 2011-2013 Mario Mulansky
+# Copyright 2012 Karsten Ahnert
+# Copyright 2013 Pascal Germroth
+# 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)
+
+use-project /boost : $(BOOST_ROOT) ;
+import openmp : * ;
+
+project
+    : requirements
+      <include>..
+      <define>BOOST_ALL_NO_LIB=1
+      <library>/boost//timer
+      [ openmp ]
+    ;
+
+exe lorenz_ensemble : lorenz_ensemble.cpp ;
+exe lorenz_ensemble_simple : lorenz_ensemble_simple.cpp ;
+exe lorenz_ensemble_nested : lorenz_ensemble_nested.cpp ;
+exe phase_chain : phase_chain.cpp ;
+exe phase_chain_omp_state : phase_chain_omp_state.cpp ;
diff --git a/examples/openmp/lorenz_ensemble.cpp b/examples/openmp/lorenz_ensemble.cpp
new file mode 100644
index 0000000..6717a50
--- /dev/null
+++ b/examples/openmp/lorenz_ensemble.cpp
@@ -0,0 +1,91 @@
+/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble.cpp
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ Parallelized Lorenz ensembles
+
+ 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 <omp.h>
+#include <vector>
+#include <iostream>
+#include <iterator>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+#include <boost/lexical_cast.hpp>
+#include "point_type.hpp"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef point<double, 3> point_type;
+typedef vector< point_type > inner_state_type;
+typedef openmp_state<point_type> state_type;
+
+const double sigma = 10.0;
+const double b = 8.0 / 3.0;
+
+
+struct sys_func {
+    const vector<double> &R;
+    sys_func( vector<double> &R ) : R(R) {}
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const {
+#       pragma omp parallel for
+        for(size_t j = 0 ; j < x.size() ; j++) {
+            size_t offset = 0;
+            for(size_t i = 0 ; i < j ; i++)
+                offset += x[i].size();
+
+            for(size_t i = 0 ; i < x[j].size() ; i++) {
+                const point_type &xi = x[j][i];
+                point_type &dxdti = dxdt[j][i];
+                dxdti[0] = -sigma * (xi[0] - xi[1]);
+                dxdti[1] = R[offset + i] * xi[0] - xi[1] - xi[0] * xi[2];
+                dxdti[2] = -b * xi[2] + xi[0] * xi[1];
+            }
+        }
+    }
+};
+
+
+int main(int argc, char **argv) {
+    size_t n = 1024;
+    if(argc > 1) n = boost::lexical_cast<size_t>(argv[1]);
+
+    vector<double> R(n);
+    const double Rmin = 0.1, Rmax = 50.0;
+#   pragma omp parallel for
+    for(size_t i = 0 ; i < n ; i++)
+        R[i] = Rmin + (Rmax - Rmin) / (n - 1) * i;
+
+    vector<point_type> inner(n, point_type(10, 10, 10));
+    state_type state;
+    split(inner, state);
+
+    cerr << "openmp_state split " << n << " into";
+    for(size_t i = 0 ; i != state.size() ; i++)
+        cerr << ' ' << state[i].size();
+    cerr << endl;
+
+    typedef runge_kutta4< state_type, double > stepper;
+
+    const double t_max = 10.0, dt = 0.01;
+
+    integrate_const(
+        stepper(),
+        sys_func(R),
+        state,
+        0.0, t_max, dt
+    );
+
+    unsplit(state, inner);
+    std::copy( inner.begin(), inner.end(), ostream_iterator<point_type>(cout, "\n") );
+
+    return 0;
+}
diff --git a/examples/openmp/lorenz_ensemble_nested.cpp b/examples/openmp/lorenz_ensemble_nested.cpp
new file mode 100644
index 0000000..4609c47
--- /dev/null
+++ b/examples/openmp/lorenz_ensemble_nested.cpp
@@ -0,0 +1,75 @@
+/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble_nested.cpp
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Pascal Germroth
+ Copyright 2013 Mario Mulansky
+
+ Parallelized Lorenz ensembles using nested omp algebra
+
+ 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 <omp.h>
+#include <vector>
+#include <iostream>
+#include <iterator>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+#include <boost/lexical_cast.hpp>
+#include "point_type.hpp"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+typedef point<double, 3> point_type;
+typedef vector< point_type > state_type;
+
+const double sigma = 10.0;
+const double b = 8.0 / 3.0;
+
+
+struct sys_func {
+    const vector<double> &R;
+    sys_func( vector<double> &R ) : R(R) {}
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const {
+#       pragma omp parallel for
+        for(size_t i = 0 ; i < x.size() ; i++) {
+            dxdt[i][0] = -sigma * (x[i][0] - x[i][1]);
+            dxdt[i][1] = R[i] * x[i][0] - x[i][1] - x[i][0] * x[i][2];
+            dxdt[i][2] = -b * x[i][2] + x[i][0] * x[i][1];
+        }
+    }
+};
+
+
+int main(int argc, char **argv) {
+    size_t n = 1024;
+    if(argc > 1) n = boost::lexical_cast<size_t>(argv[1]);
+
+    vector<double> R(n);
+    const double Rmin = 0.1, Rmax = 50.0;
+#   pragma omp parallel for
+    for(size_t i = 0 ; i < n ; i++)
+        R[i] = Rmin + (Rmax - Rmin) / (n - 1) * i;
+
+    state_type state( n , point_type(10, 10, 10) );
+
+    typedef runge_kutta4< state_type, double , state_type , double ,
+                          openmp_nested_algebra<vector_space_algebra> > stepper;
+
+    const double t_max = 10.0, dt = 0.01;
+
+    integrate_const(
+        stepper(),
+        sys_func(R),
+        state,
+        0.0, t_max, dt
+    );
+
+    std::copy( state.begin(), state.end(), ostream_iterator<point_type>(cout, "\n") );
+
+    return 0;
+}
diff --git a/examples/openmp/lorenz_ensemble_simple.cpp b/examples/openmp/lorenz_ensemble_simple.cpp
new file mode 100644
index 0000000..fbaf249
--- /dev/null
+++ b/examples/openmp/lorenz_ensemble_simple.cpp
@@ -0,0 +1,79 @@
+/* Boost libs/numeric/odeint/examples/openmp/lorenz_ensemble_simple.cpp
+
+ Copyright 2013 Karsten Ahnert
+ Copyright 2013 Mario Mulansky
+ Copyright 2013 Pascal Germroth
+
+ Parallelized Lorenz ensembles
+
+ 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 <omp.h>
+#include <vector>
+#include <iostream>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+#include "point_type.hpp"
+
+using namespace std;
+
+typedef vector<double> vector_type;
+typedef point<double, 3> point_type;
+typedef vector<point_type> state_type;
+
+const double sigma = 10.0;
+const double b = 8.0 / 3.0;
+
+struct sys_func {
+    const vector_type &R;
+    sys_func( const vector_type &_R ) : R( _R ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double t ) const {
+        const size_t n = x.size();
+#       pragma omp parallel for
+        for(size_t i = 0 ; i < n ; i++) {
+            const point_type &xi = x[i];
+            point_type &dxdti = dxdt[i];
+            dxdti[0] = -sigma * (xi[0] - xi[1]);
+            dxdti[1] = R[i] * xi[0] - xi[1] - xi[0] * xi[2];
+            dxdti[2] = -b * xi[2] + xi[0] * xi[1];
+        }
+    }
+};
+
+
+int main() {
+    using namespace boost::numeric::odeint;
+
+    const size_t n = 1024;
+    vector_type R(n);
+    const double Rmin = 0.1, Rmax = 50.0;
+#   pragma omp parallel for
+    for(size_t i = 0 ; i < n ; i++)
+        R[i] = Rmin + (Rmax - Rmin) / (n - 1) * i;
+
+    state_type X(n, point_type(10, 10, 10));
+
+    const double t_max = 10.0, dt = 0.01;
+
+    // Simple stepper with constant step size
+    // typedef runge_kutta4<state_type, double, state_type, double,
+    //                      openmp_range_algebra> stepper;
+
+    // integrate_const(stepper(), sys_func(R), X, 0.0, t_max, dt);
+
+    // Controlled stepper with variable step size
+    typedef runge_kutta_fehlberg78<state_type, double, state_type, double,
+                                   openmp_range_algebra> error_stepper_type;
+    typedef controlled_runge_kutta<error_stepper_type> controlled_stepper_type;
+    controlled_stepper_type controlled_stepper;
+
+    integrate_adaptive(controlled_stepper, sys_func(R), X, 0.0, t_max, dt);
+
+    copy( X.begin(), X.end(), ostream_iterator<point_type>(cout, "\n") );
+
+    return 0;
+}
diff --git a/examples/openmp/openmp.jam b/examples/openmp/openmp.jam
new file mode 100644
index 0000000..226bcf5
--- /dev/null
+++ b/examples/openmp/openmp.jam
@@ -0,0 +1,52 @@
+# Copyright 2013 Karsten Ahnert
+# Copyright 2013 Mario Mulansky
+# Copyright 2013 Pascal Germroth
+# 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)
+
+
+# Only builds target with supported OpenMP enabled toolsets.
+#
+# use as:
+#   exe omp : omp.cpp : [ openmp ] ;
+#
+rule openmp return
+    # default
+    <build>no
+    # GNU C++
+    <toolset>gcc:<cxxflags>-fopenmp
+    <toolset>gcc:<linkflags>-fopenmp
+    <toolset>gcc:<build>yes
+    # Microsoft Visual C++
+    <toolset>msvc:<cxxflags>/openmp
+    <toolset>msvc:<linkflags>/openmp
+    <toolset>msvc:<build>yes
+    # Intel C++
+    <toolset>intel-linux:<cxxflags>-openmp
+    <toolset>intel-linux:<linkflags>-openmp
+    <toolset>intel-linux:<build>yes
+    <toolset>intel-win:<cxxflags>-Qopenmp
+    <toolset>intel-win:<linkflags>-Qopenmp
+    <toolset>intel-win:<build>yes
+    # HP aC++
+    <toolset>acc:<cxxflags>+Oopenmp
+    <toolset>acc:<linkflags>+Oopenmp
+    <toolset>acc:<build>yes
+    # Sun Studio
+    <toolset>sun:<cxxflags>-xopenmp
+    <toolset>sun:<linkflags>-xopenmp
+    <toolset>sun:<build>yes
+    # IBM XL
+    <toolset>vacpp:<cxxflags>-qsmp=omp
+    <toolset>vacpp:<linkflags>-qsmp=omp
+    <toolset>vacpp:<build>yes
+    # PG++
+    <toolset>pgi:<cxxflags>-mp
+    <toolset>pgi:<linkflags>-mp
+    <toolset>pgi:<build>yes
+    # Pathscale
+    <toolset>pathscale:<cxxflags>-mp
+    <toolset>pathscale:<linkflags>-mp
+    <toolset>pathscale:<build>yes
+    ;
diff --git a/examples/openmp/phase_chain.cpp b/examples/openmp/phase_chain.cpp
new file mode 100644
index 0000000..08876fc
--- /dev/null
+++ b/examples/openmp/phase_chain.cpp
@@ -0,0 +1,95 @@
+/*
+ * phase_chain.cpp
+ *
+ * Example of OMP parallelization with odeint
+ *
+ * Copyright 2013 Karsten Ahnert
+ * Copyright 2013 Mario Mulansky
+ * Copyright 2013 Pascal Germroth
+ * 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 <vector>
+#include <boost/random.hpp>
+#include <boost/timer/timer.hpp>
+//[phase_chain_openmp_header
+#include <omp.h>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+//]
+
+using namespace std;
+using namespace boost::numeric::odeint;
+using boost::timer::cpu_timer;
+using boost::math::double_constants::pi;
+
+//[phase_chain_vector_state
+typedef std::vector< double > state_type;
+//]
+
+//[phase_chain_rhs
+struct phase_chain
+{
+    phase_chain( double gamma = 0.5 )
+    : m_gamma( gamma ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
+    {
+        const size_t N = x.size();
+        #pragma omp parallel for schedule(runtime)
+        for(size_t i = 1 ; i < N - 1 ; ++i)
+        {
+            dxdt[i] = coupling_func( x[i+1] - x[i] ) +
+                      coupling_func( x[i-1] - x[i] );
+        }
+        dxdt[0  ] = coupling_func( x[1  ] - x[0  ] );
+        dxdt[N-1] = coupling_func( x[N-2] - x[N-1] );
+    }
+
+    double coupling_func( double x ) const
+    {
+        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
+    }
+
+    double m_gamma;
+};
+//]
+
+
+int main( int argc , char **argv )
+{
+    //[phase_chain_init
+    size_t N = 131101;
+    state_type x( N );
+    boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
+    boost::random::mt19937 engine( 0 );
+    generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
+    //]
+
+    //[phase_chain_stepper
+    typedef runge_kutta4<
+                      state_type , double ,
+                      state_type , double ,
+                      openmp_range_algebra
+                    > stepper_type;
+    //]
+
+    //[phase_chain_scheduling
+    int chunk_size = N/omp_get_max_threads();
+    omp_set_schedule( omp_sched_static , chunk_size );
+    //]
+
+    cpu_timer timer;
+    //[phase_chain_integrate
+    integrate_n_steps( stepper_type() , phase_chain( 1.2 ) ,
+                       x , 0.0 , 0.01 , 100 );
+    //]
+    double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
+    std::cerr << run_time << "s" << std::endl;
+    // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
+
+    return 0;
+}
diff --git a/examples/openmp/phase_chain_omp_state.cpp b/examples/openmp/phase_chain_omp_state.cpp
new file mode 100644
index 0000000..8627f58
--- /dev/null
+++ b/examples/openmp/phase_chain_omp_state.cpp
@@ -0,0 +1,98 @@
+/*
+ * phase_chain_omp_state.cpp
+ *
+ * Example of OMP parallelization with odeint
+ *
+ * Copyright 2013 Karsten Ahnert
+ * Copyright 2013 Mario Mulansky
+ * Copyright 2013 Pascal Germroth
+ * 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 <vector>
+#include <boost/random.hpp>
+#include <boost/timer/timer.hpp>
+#include <omp.h>
+#include <boost/numeric/odeint.hpp>
+#include <boost/numeric/odeint/external/openmp/openmp.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+using namespace boost::numeric::odeint;
+using boost::timer::cpu_timer;
+using boost::math::double_constants::pi;
+
+typedef openmp_state<double> state_type;
+
+//[phase_chain_state_rhs
+struct phase_chain_omp_state
+{
+    phase_chain_omp_state( double gamma = 0.5 )
+    : m_gamma( gamma ) { }
+
+    void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const
+    {
+        const size_t N = x.size();
+        #pragma omp parallel for schedule(runtime)
+        for(size_t n = 0 ; n < N ; ++n)
+        {
+            const size_t M = x[n].size();
+            for(size_t m = 1 ; m < M-1 ; ++m)
+            {
+                dxdt[n][m] = coupling_func( x[n][m+1] - x[n][m] ) +
+                             coupling_func( x[n][m-1] - x[n][m] );
+            }
+            dxdt[n][0] = coupling_func( x[n][1] - x[n][0] );
+            if( n > 0 )
+            {
+                dxdt[n][0] += coupling_func( x[n-1].back() - x[n].front() );
+            }
+            dxdt[n][M-1] = coupling_func( x[n][M-2] - x[n][M-1] );
+            if( n < N-1 )
+            {
+                dxdt[n][M-1] += coupling_func( x[n+1].front() - x[n].back() );
+            }
+        }
+    }
+
+    double coupling_func( double x ) const
+    {
+        return sin( x ) - m_gamma * ( 1.0 - cos( x ) );
+    }
+
+    double m_gamma;
+};
+//]
+
+
+int main( int argc , char **argv )
+{
+    //[phase_chain_state_init
+    const size_t N = 131101;
+    vector<double> x( N );
+    boost::random::uniform_real_distribution<double> distribution( 0.0 , 2.0*pi );
+    boost::random::mt19937 engine( 0 );
+    generate( x.begin() , x.end() , boost::bind( distribution , engine ) );
+    const size_t blocks = omp_get_max_threads();
+    state_type x_split( blocks );
+    split( x , x_split );
+    //]
+
+
+    cpu_timer timer;
+    //[phase_chain_state_integrate
+    integrate_n_steps( runge_kutta4<state_type>() , phase_chain_omp_state( 1.2 ) ,
+                       x_split , 0.0 , 0.01 , 100 );
+    unsplit( x_split , x );
+    //]
+
+    double run_time = static_cast<double>(timer.elapsed().wall) * 1.0e-9;
+    std::cerr << run_time << "s" << std::endl;
+    // copy(x.begin(), x.end(), ostream_iterator<double>(cout, "\n"));
+
+    return 0;
+}