diff --git a/examples/2d_lattice/spreading.cpp b/examples/2d_lattice/spreading.cpp
new file mode 100644
index 0000000..88c60bf
--- /dev/null
+++ b/examples/2d_lattice/spreading.cpp
@@ -0,0 +1,122 @@
+/*
+ Copyright 2011 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)
+ */
+
+
+/*
+ * Example of a 2D simulation of nonlinearly coupled oscillators.
+ * Program just prints final energy which should be close to the initial energy (1.0).
+ * No parallelization is employed here.
+ * Run time on a 2.3GHz Intel Core-i5: about 10 seconds for 100 steps.
+ * Compile simply via bjam or directly:
+ * g++ -O3 -I${BOOST_ROOT} -I../../../../.. spreading.cpp
+ */
+
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cstdlib>
+#include <sys/time.h>
+
+#include <boost/ref.hpp>
+#include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp>
+
+// we use a vector< vector< double > > as state type,
+// for that some functionality has to be added for odeint to work
+#include "nested_range_algebra.hpp"
+#include "vector_vector_resize.hpp"
+
+// defines the rhs of our dynamical equation
+#include "lattice2d.hpp"
+/* dynamical equations (Hamiltonian structure):
+dqdt_{i,j} = p_{i,j}
+dpdt_{i,j} = - omega_{i,j}*q_{i,j} - \beta*[ (q_{i,j} - q_{i,j-1})^3
+                                            +(q_{i,j} - q_{i,j+1})^3
+                                            +(q_{i,j} - q_{i-1,j})^3
+                                            +(q_{i,j} - q_{i+1,j})^3 ]
+*/
+
+
+using namespace std;
+
+static const int MAX_N = 1024;//2048;
+
+static const size_t KAPPA = 2;
+static const size_t LAMBDA = 4;
+static const double W = 1.0;
+static const double gap = 0.0;
+static const size_t steps = 100;
+static const double dt = 0.1;
+
+double initial_e = 1.0;
+double beta = 1.0;
+int realization_index = 0;
+
+//the state type
+typedef vector< vector< double > > state_type;
+
+//the stepper, choose a symplectic one to account for hamiltonian structure
+//use nested_range_algebra for calculations on vector< vector< ... > >
+typedef boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan< 
+    state_type , state_type , double , state_type , state_type , double , 
+    nested_range_algebra< boost::numeric::odeint::range_algebra > ,
+    boost::numeric::odeint::default_operations > stepper_type;
+
+double time_diff_in_ms( timeval &t1 , timeval &t2 )
+{ return (t2.tv_sec - t1.tv_sec)*1000.0 + (t2.tv_usec - t1.tv_usec)/1000.0 + 0.5; }
+
+
+int main( int argc, const char* argv[] ) {
+
+    srand( time(NULL) );
+
+    lattice2d< KAPPA , LAMBDA > lattice( beta );
+
+
+    lattice.generate_pot( W , gap , MAX_N );
+
+    state_type q( MAX_N , vector< double >( MAX_N , 0.0 ) );
+
+    state_type p( q );
+
+    state_type energy( q );
+
+    p[MAX_N/2][MAX_N/2] = sqrt( 0.5*initial_e );
+    p[MAX_N/2+1][MAX_N/2] = sqrt( 0.5*initial_e );
+    p[MAX_N/2][MAX_N/2+1] = sqrt( 0.5*initial_e );
+    p[MAX_N/2+1][MAX_N/2+1] = sqrt( 0.5*initial_e );
+
+    cout.precision(10);
+
+    lattice.local_energy( q , p , energy );
+    double e=0.0;
+    for( size_t i=0 ; i<energy.size() ; ++i )
+        for( size_t j=0 ; j<energy[i].size() ; ++j )
+        {
+            e += energy[i][j];
+        }
+
+    cout << "initial energy: " << lattice.energy( q , p ) << endl;
+
+    timeval elapsed_time_start , elapsed_time_end;
+    gettimeofday(&elapsed_time_start , NULL);
+
+    stepper_type stepper;
+
+    for( size_t step=0 ; step<=steps ; ++step )
+    {
+        stepper.do_step( boost::ref( lattice ) , 
+                         make_pair( boost::ref( q ) , boost::ref( p ) ) , 
+                         0.0 , 0.1 );
+    }
+
+    gettimeofday(&elapsed_time_end , NULL);
+    double elapsed_time = time_diff_in_ms( elapsed_time_start , elapsed_time_end );
+    cout << steps << " steps in " << elapsed_time/1000 << " s (energy: " << lattice.energy( q , p ) << ")" << endl;
+}