diff --git a/doc/concepts/controlled_stepper.qbk b/doc/concepts/controlled_stepper.qbk
new file mode 100644
index 0000000..20afa85
--- /dev/null
+++ b/doc/concepts/controlled_stepper.qbk
@@ -0,0 +1,74 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright 2011-2012 Karsten Ahnert
+  Copyright 2011 Mario Mulansky
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+[section Controlled Stepper]
+
+This concept specifies the interface a controlled stepper has to fulfill to be used within __integrate_functions.
+
+[heading Description]
+
+A controlled stepper following this Controlled Stepper concept provides the possibility to perform one step of the solution /x(t)/ of an ODE with step-size /dt/ to obtain /x(t+dt)/ with a given step-size /dt/.
+Depending on an error estimate of the solution the step might be rejected and a smaller step-size is suggested.
+
+[heading Associated types]
+
+*  '''<para>'''[*state_type]'''</para>'''
+'''<para>'''`Stepper::state_type`'''</para>'''
+'''<para>'''The type characterizing the state of the ODE, hence ['x].'''</para>'''
+
+*  '''<para>'''[*deriv_type]'''</para>'''
+'''<para>'''`Stepper::deriv_type`'''</para>'''
+'''<para>'''The type characterizing the derivative of the ODE, hence ['d x/dt].'''</para>'''
+
+*  '''<para>'''[*time_type]'''</para>'''
+'''<para>'''`Stepper::time_type`'''</para>'''
+'''<para>'''The type characterizing the dependent variable of the ODE, hence the time ['t].'''</para>'''
+
+*  '''<para>'''[*value_type]'''</para>'''
+'''<para>'''`Stepper::value_type`'''</para>'''
+'''<para>'''The numerical data type which is used within the stepper, something like `float`, `double`, `complex&lt; double &gt;`.'''</para>'''
+
+*  '''<para>'''[*stepper_category]'''</para>'''
+'''<para>'''`Stepper::stepper_category`'''</para>'''
+'''<para>'''A tag type characterizing the category of the stepper. This type must be convertible to `controlled_stepper_tag`.'''</para>'''
+
+
+
+[heading Notation]
+
+[variablelist 
+  [[`ControlledStepper`] [A type that is a model of Controlled Stepper]]
+  [[`State`] [A type representing the state /x/ of the ODE]]
+  [[`Time`] [A type representing the time /t/ of the ODE]]
+  [[`stepper`] [An object of type `ControlledStepper`]]
+  [[`x`] [Object of type `State`]]
+  [[`t`, `dt`] [Objects of type `Time`]]
+  [[`sys`] [An object defining the ODE, should be a model of __system, __symplectic_system, __simple_symplectic_system or __implicit_system.]]
+]
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Do step] [``stepper.try_step( sys , x , t , dt )``] [`controlled_step_result`] [Tries one step of step size `dt`. If the step was successful, `success` is returned, the resulting state is written to `x`, the new time is stored in `t` and `dt` now contains a new (possibly larger) step-size for the next step. If the error was too big, `rejected` is returned and the results are neglected - `x` and `t` are unchanged and `dt` now contains a reduced step-size to be used for the next try.] ]
+  [/ [Do step with reference] [`stepper.try_step( boost::ref(sys) , x , t , dt )`] [`void`] [Same as above with `System` as reference] ]
+]
+
+[heading Models]
+
+* `controlled_error_stepper< runge_kutta_cash_karp54 >`
+* `controlled_error_stepper_fsal< runge_kutta_dopri5 >`
+* `controlled_error_stepper< runge_kutta_fehlberg78 >`
+* `rosenbrock4_controller`
+* `bulirsch_stoer`
+
+[endsect]
\ No newline at end of file
diff --git a/doc/concepts/dense_output_stepper.qbk b/doc/concepts/dense_output_stepper.qbk
new file mode 100644
index 0000000..19d3a25
--- /dev/null
+++ b/doc/concepts/dense_output_stepper.qbk
@@ -0,0 +1,85 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright (c) 2009-2013 Karsten Ahnert
+  Copyright (c) 2009-2013 Mario Mulansky
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+
+[section Dense Output Stepper]
+
+This concept specifies the interface a dense output stepper has to fulfill to be used  within __integrate_functions.
+
+[heading Description]
+A dense output stepper following this Dense Output Stepper concept provides the possibility to perform a single step of the solution /x(t)/ of an ODE to obtain /x(t+dt)/.
+The step-size `dt` might be adjusted automatically due to error control.
+Dense output steppers also can interpolate the solution to calculate the state /x(t')/ at any point /t <= t' <= t+dt/.
+
+[heading Associated types]
+
+*  '''<para>'''[*state_type]'''</para>'''
+'''<para>'''`Stepper::state_type`'''</para>'''
+'''<para>'''The type characterizing the state of the ODE, hence ['x].'''</para>'''
+
+*  '''<para>'''[*deriv_type]'''</para>'''
+'''<para>'''`Stepper::deriv_type`'''</para>'''
+'''<para>'''The type characterizing the derivative of the ODE, hence ['d x/dt].'''</para>'''
+
+*  '''<para>'''[*time_type]'''</para>'''
+'''<para>'''`Stepper::time_type`'''</para>'''
+'''<para>'''The type characterizing the dependent variable of the ODE, hence the time ['t].'''</para>'''
+
+*  '''<para>'''[*value_type]'''</para>'''
+'''<para>'''`Stepper::value_type`'''</para>'''
+'''<para>'''The numerical data type which is used within the stepper, something like `float`, `double`, `complex&lt; double &gt;`.'''</para>'''
+
+*  '''<para>'''[*stepper_category]'''</para>'''
+'''<para>'''`Stepper::stepper_category`'''</para>'''
+'''<para>'''A tag type characterizing the category of the stepper. This type must be convertible to `dense_output_stepper_tag`.'''</para>'''
+
+
+[heading Notation]
+
+[variablelist 
+  [[`Stepper`] [A type that is a model of Dense Output Stepper]]
+  [[`State`] [A type representing the state /x/ of the ODE]]
+  [[`stepper`] [An object of type `Stepper`]]
+  [[`x0`, `x`] [Object of type `State`]]
+  [[`t0`, `dt0`, `t`] [Objects of type `Stepper::time_type`]]
+  [[`sys`] [An object defining the ODE, should be a model of __system, __symplectic_system, __simple_symplectic_system or __implicit_system.]]
+]
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+
+  [[Initialize integration] [`stepper.initialize( x0 , t0 , dt0 )`] [void] [Initializes the stepper with initial values `x0`, `t0` and `dt0`.]]
+
+  [[Do step] [`stepper.do_step( sys )`] [`std::pair< Stepper::time_type , Stepper::time_type >`] [Performs one step using the ODE defined by `sys`. The step-size might be changed internally due to error control. This function returns a pair containing `t` and `t+dt` representing the interval for which interpolation can be performed.] ]
+
+  [/ [Do step with reference] [`stepper.do_step( boost::ref( sys ) )`] [`std::pair< Stepper::time_type , Stepper::time_type >`] [Same as above with `System` as reference] ]
+
+  [[Do interpolation] [`stepper.calc_state( t_inter , x )`] [`void`] [Performs the interpolation to calculate /x(t[sub inter]/) where /t <= t[sub inter] <= t+dt/.]]
+
+  [[Get current time] [`stepper.current_time()`] [`const Stepper::time_type&`] [Returns the current time /t+dt/ of the stepper, that is the end time of the last step and the starting time for the next call of `do_step`]]
+
+  [[Get current state] [`stepper.current_state()`] [`const Stepper::state_type&`] [Returns the current state of the stepper, that is /x(t+dt)/, the state at the time returned by `stepper.current_time()`]] 
+
+  [[Get current time step] [`stepper.current_time_step()`] [`const
+  Stepper::time_type&`] [Returns the current step size of the stepper, that is
+  /dt/]] 
+]
+
+[heading Models]
+
+* `dense_output_controlled_explicit_fsal< controlled_error_stepper_fsal< runge_kutta_dopri5 >`
+* `bulirsch_stoer_dense_out`
+* `rosenbrock4_dense_output`
+
+[endsect]
diff --git a/doc/concepts/error_stepper.qbk b/doc/concepts/error_stepper.qbk
new file mode 100644
index 0000000..8ebfde9
--- /dev/null
+++ b/doc/concepts/error_stepper.qbk
@@ -0,0 +1,101 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright 2011-2013 Karsten Ahnert
+  Copyright 2011 Mario Mulansky
+  Copyright 2012 Sylwester Arabas
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+
+[section Error Stepper]
+
+This concepts specifies the interface an error stepper has to fulfill to be used within a ControlledErrorStepper. An error stepper must always fulfill the stepper concept. This can trivially implemented by
+
+``
+template< class System >
+error_stepper::do_step( System sys , state_type &x , time_type t , time_type dt )
+{
+    state_type xerr;
+    // allocate xerr
+    do_step( sys , x , t , dt , xerr );
+}
+``
+
+[heading Description]
+
+An error stepper following this Error Stepper concept is capable of doing one step of the solution /x(t)/ of an ODE with step-size /dt/ to obtain /x(t+dt)/ and also computing an error estimate ['x[sub err]] of the result.
+Error Steppers can be Runge-Kutta steppers, symplectic steppers as well as implicit steppers.
+Based on the stepper type, the ODE is defined as __system, __symplectic_system, __simple_symplectic_system or __implicit_system.
+
+[heading Refinement of]
+
+* DefaultConstructable
+* CopyConstructable
+* Stepper
+
+[heading Associated types]
+
+*  '''<para>'''[*state_type]'''</para>'''
+'''<para>'''`Stepper::state_type`'''</para>'''
+'''<para>'''The type characterizing the state of the ODE, hence ['x].'''</para>'''
+
+*  '''<para>'''[*deriv_type]'''</para>'''
+'''<para>'''`Stepper::deriv_type`'''</para>'''
+'''<para>'''The type characterizing the derivative of the ODE, hence ['d x/dt].'''</para>'''
+
+*  '''<para>'''[*time_type]'''</para>'''
+'''<para>'''`Stepper::time_type`'''</para>'''
+'''<para>'''The type characterizing the dependent variable of the ODE, hence the time ['t].'''</para>'''
+
+*  '''<para>'''[*value_type]'''</para>'''
+'''<para>'''`Stepper::value_type`'''</para>'''
+'''<para>'''The numerical data type which is used within the stepper, something like `float`, `double`, `complex&lt; double &gt;`.'''</para>'''
+
+*  '''<para>'''[*order_type]'''</para>'''
+'''<para>'''`Stepper::order_type`'''</para>'''
+'''<para>'''The type characterizing the order of the ODE, typically `unsigned short`.'''</para>'''
+
+*  '''<para>'''[*stepper_category]'''</para>'''
+'''<para>'''`Stepper::stepper_category`'''</para>'''
+'''<para>'''A tag type characterizing the category of the stepper. This type must be convertible to `error_stepper_tag`.'''</para>'''
+
+
+[heading Notation]
+
+[variablelist 
+  [[`ErrorStepper`] [A type that is a model of Error Stepper]]
+  [[`State`] [A type representing the state /x/ of the ODE]]
+  [[`Error`] [A type representing the error calculated by the stepper, usually same as `State`]]
+  [[`Time`] [A type representing the time /t/ of the ODE]]
+  [[`stepper`] [An object of type `ErrorStepper`]]
+  [[`x`] [Object of type `State`]]
+  [[`xerr`] [Object of type `Error`]]
+  [[`t`, `dt`] [Objects of type `Time`]]
+  [[`sys`] [An object defining the ODE, should be a model of either __system, __symplectic_system, __simple_symplectic_system or __implicit_system.]]
+]
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Get the stepper order] [`stepper.order()`] [`order_type`] [Returns the order of the stepper for one step without error estimation.]]
+  [[Get the stepper order] [`stepper.stepper_order()`] [`order_type`] [Returns the order of the stepper for one error estimation step which is used for error calculation.]]
+  [[Get the error order] [`stepper.errorr_order()`] [`order_type`] [Returns the order of the error step which is used for error calculation.]]
+  [[Do step] [`stepper.do_step( sys , x , t , dt )`] [`void`] [Performs one step of step size `dt`. The newly obtained state is written in-place to `x`.] ]
+  [[Do step with error estimation] [`stepper.do_step( sys , x , t , dt , xerr )`] [`void`] [Performs one step of step size `dt` with error estimation.  The newly obtained state is written in-place to `x` and the estimated error to `xerr`.] ]
+  [/ [Do step with reference] [`stepper.do_step( boost::ref(sys) , x , t , dt , xerr )`] [`void`] [Performs one step of step size `dt`.  The newly obtained state is written in-place to `x` and the estimated error to `xerr`.] ]
+]
+
+[heading Models]
+
+* `runge_kutta_cash_karp54`
+* `runge_kutta_dopri5`
+* `runge_kutta_fehlberg78`
+* `rosenbrock4`
+
+[endsect]
diff --git a/doc/concepts/implicit_system.qbk b/doc/concepts/implicit_system.qbk
new file mode 100644
index 0000000..e889023
--- /dev/null
+++ b/doc/concepts/implicit_system.qbk
@@ -0,0 +1,43 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright 2011 Mario Mulansky
+  Copyright 2011-2012 Karsten Ahnert
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+
+[section Implicit System]
+
+[heading Description]
+
+This concept describes how to define a ODE that can be solved by an implicit routine.
+Implicit routines need not only the function /f(x,t)/ but also the Jacobian /df/dx = A(x,t)/.
+/A/ is a matrix and implicit routines need to solve the linear problem /Ax = b/.
+In odeint this is implemented with use of __ublas, therefore, the ['state_type] implicit routines is ['ublas::vector] and the matrix is defined as ['ublas::matrix].
+
+[heading Notation]
+
+[variablelist
+  [[`System`] [A type that is a model of `Implicit System`]]
+  [[`Time`] [A type representing the time of the ODE]]
+  [[`sys`] [An object of type `System`]]
+  [[`x`] [Object of type ublas::vector]]
+  [[`dxdt`] [Object of type ublas::vector]]
+  [[`jacobi`] [Object of type ublas::matrix]]
+  [[`t`] [Object of type `Time`]]
+]
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Calculate ['dx/dt := f(x,t)]] [`sys.first( x , dxdt , t )`] [`void`] [Calculates `f(x,t)`, the result is stored into dxdt] ]
+  [[Calculate ['A := df/dx (x,t)]] [`sys.second( x , jacobi , t )`] [`void`] [Calculates the Jacobian of /f/ at /x/,/t/, the result is stored into `jacobi`] ]
+]
+
+[endsect]
\ No newline at end of file
diff --git a/doc/concepts/second_order_system.qbk b/doc/concepts/second_order_system.qbk
new file mode 100644
index 0000000..9ca9009
--- /dev/null
+++ b/doc/concepts/second_order_system.qbk
@@ -0,0 +1,44 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright (c) 2009-2013 Karsten Ahnert
+  Copyright (c) 2009-2013 Mario Mulansky
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+[section Second Order System]
+
+[heading Description]
+
+The Second Order System concept models the algorithmic implementation of the rhs for steppers requirering the second order
+derivative, hence the r.h.s. of the ODE ['x'' = f(x,x',t)]. The only requirement for this concept is that it should be callable
+with a specific parameter syntax (see below). A Second Order System is typically implemented as a function or a functor.
+Systems fulfilling this concept are required by the Velocity Verlet method.
+
+[heading Notation]
+
+[variablelist
+  [[`System`] [A type that is a model of Second Order System]] 
+  [[`Space`] [A type representing the state /x/ of the ODE]]
+  [[`Velocity`] [A type representing the derivative /x'/ of the ODE]]
+  [[`Acceleration`] [A type representing the second order derivative /x''/ of the ODE]]
+  [[`Time`] [A type representing the time]]
+  [[`sys`] [An object of type `System`]]
+  [[`x`] [Object of type `Space`]]
+  [[`v`] [Object of type `Velocity`]]
+  [[`a`] [Object of type `Acceleration`]]
+  [[`t`] [Object of type `Time`]]
+]
+
+[heading Valid expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Calculate ['x'' := f(x,x',t)]] [`sys( x , v , a , t )`] [`void`] [Calculates f(x,x',t), the result is stored into a.] ]
+]
+
+[endsect]
\ No newline at end of file
diff --git a/doc/concepts/state_algebra_operations.qbk b/doc/concepts/state_algebra_operations.qbk
new file mode 100644
index 0000000..d92bb83
--- /dev/null
+++ b/doc/concepts/state_algebra_operations.qbk
@@ -0,0 +1,126 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright 2011 Mario Mulansky
+  Copyright 2012 Karsten Ahnert
+  Copyright 2013 Pascal Germroth
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+
+[section State Algebra Operations]
+
+[note The following does not apply to implicit steppers like implicit_euler or Rosenbrock 4 as there the `state_type` can not be changed from `ublas::vector` and no algebra/operations are used.]
+
+[heading Description]
+
+The `State`, `Algebra` and `Operations` together define a concept describing how the mathematical vector operations required for the stepper algorithms are performed.
+The typical vector operation done within steppers is 
+
+['*y* = __Sigma __alpha[sub i] [*x[sub i]]].
+
+The `State` represents the state variable of an ODE, usually denoted with /x/.
+Algorithmically, the state is often realized as a `vector< double >` or `array< double , N >`, however, the genericity of odeint enables you to basically use anything as a state type.
+The algorithmic counterpart of such mathematical expressions is divided into two parts.
+First, the `Algebra` is used to account for the vector character of the equation.
+In the case of a `vector` as state type this means the `Algebra` is responsible for iteration over all vector elements.
+Second, the `Operations` are used to represent the actual operation applied to each of the vector elements.
+So the `Algebra` iterates over all elements of the `State`s and calls an operation taken from the `Operations` for each element.
+This is where `State`, `Algebra` and `Operations` have to work together to make odeint running.
+Please have a look at the `range_algebra` and `default_operations` to see an example how this is implemented.
+
+In the following we describe how `State`, `Algebra` and `Operations` are used together within the stepper implementations.
+
+[section Operations]
+
+[heading Notation]
+
+[variablelist
+  [[`Operations`] [The operations type]]
+  [/[`Time`] [A type representing the time type of steppers]]
+  [[`Value1`, ... , `ValueN`] [Types representing the value or time type of stepper]]
+  [[`Scale`] [Type of the scale operation]]
+  [[`scale`] [Object of type `Scale`]]
+  [[[^ScaleSum['N]]] [Type that represents a general scale_sum operation, [^/N/] should be replaced by a number from 1 to 14.]]
+  [[[^scale_sum['N]]] [Object of type [^ScaleSum['N]], [^/N/] should be replaced by a number from 1 to 14.]]
+  [[`ScaleSumSwap2`] [Type of the scale sum swap operation]]
+  [[`scale_sum_swap2`] [Object of type `ScaleSumSwap2`]]
+  [[`a1, a2, ...`] [Objects of type `Value1`, `Value2`, ...]]
+  [[`y, x1, x2, ...`] [Objects of `State`'s value type]]
+]
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Get scale operation] [`Operations::scale< Value >`] [`Scale`] [Get `Scale` from `Operations`]]
+  [[`Scale` constructor] [`Scale< Value >( a )`] [`Scale`] [Constructs a `Scale` object]]
+  [[`Scale` operation] [`scale( x )`] [`void`] [Calculates `x *= a`]]
+  [[Get general `scale_sum` operation] [[^Operations::scale_sum['N]< Value1 , ... , ValueN >]] [[^ScaleSum['N]]] [Get the [^ScaleSum['N]] type from `Operations`, [^/N/] should be replaced by a number from 1 to 14.]]
+  [[`scale_sum` constructor] [[^ScaleSum['N]< Value1 , ... , ValueN >( a1 , ... , aN )]] [[^ScaleSum['N]]] [Constructs a `scale_sum` object given [^/N/] parameter values with [^/N/] between 1 and 14.]]
+  [[`scale_sum` operation] [[^scale_sum['N]( y , x1 , ... , xN )]] [`void`] [Calculates `y = a1*x1 + a2*x2 + ... + aN*xN`. Note that this is an [^/N/+1]-ary function call.]]
+  [[Get scale sum swap operation] [`Operations::scale_sum_swap2< Value1 , Value2 >`] [`ScaleSumSwap2`] [Get scale sum swap from operations]]
+  [[`ScaleSumSwap2` constructor] [`ScaleSumSwap2< Value1 , Value2 >( a1 , a2 )`] [`ScaleSumSwap2`] [Constructor]]
+  [[`ScaleSumSwap2` operation] [`scale_sum_swap2( x1 , x2 , x3 )`] [`void`] [Calculates `tmp = x1`, `x1 = a1*x2 + a2*x3` and `x2 = tmp`.]]
+]
+
+[endsect]
+
+[section Algebra]
+
+[heading Notation]
+
+[variablelist
+  [[`State`] [The state type]]
+  [[`Algebra`] [The algebra type]]
+  [[[^Operation['N]]] [An [^/N/]-ary operation type, [^/N/] should be a number from 1 to 14.]]
+  [[`algebra`] [Object of type `Algebra`]]
+  [[[^operation['N]]] [Object of type [^Operation['N]]]]
+  [[`y, x1, x2, ...`] [Objects of type `State`]]
+]
+
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Vector Operation with arity 2] [`algebra.for_each2( y , x , operation2 )`] [void] [Calls `operation2( y_i , x_i )` for each element `y_i` of `y` and `x_i` of `x`.]]
+  [[Vector Operation with arity 3] [`algebra.for_each3( y , x1 , x2 , operation3 )`] [void] [Calls `operation3( y_i , x1_i , x2_i )` for each element `y_i` of `y` and `x1_i` of `x1` and `x2_i` of `x2`.]]
+  [[Vector Operation with arity [^/N/]] [[^algebra.for_each['N]( y , x1 , ... , xN , operation['N] )]] [void] [Calls [^operation['N]( y_i , x1_i , ... , xN_i )] for each element `y_i` of `y` and `x1_i` of `x1` and so on. [^/N/] should be replaced by a number between 1 and 14.]]
+]
+
+[endsect]
+
+[section Pre-Defined implementations]
+
+As standard configuration odeint uses the `range_algebra` and `default_operations` which suffices most situations.
+However, a few more possibilities exist either to gain better performance or to ensure interoperability with other libraries.
+In the following we list the existing `Algebra`/`Operations` configurations that can be used in the steppers.
+
+[table
+  [[`State`] [`Algebra`] [`Operations`] [Remarks]]
+  [[Anything supporting __boost_range, like `std::vector`, `std::list`, `boost::array`,... based on a `value_type` that supports operators +,* (typically `double`)] [`range_algebra`] [`default_operations`] [Standard implementation, applicable for most typical situations.]]
+  [[`boost::array` based on a `value_type` that supports operators +,*] [`array_algebra`] [`default_operations`] [Special implementation for boost::array with better performance than `range_algebra`]]
+  [[Anything that defines operators + within itself and * with scalar (Mathematically spoken, anything that is a vector space).] [`vector_space_algebra`] [`default_operations`] [For the use of __controlled_stepper, the template `vector_space_reduce` has to be instantiated.]]
+  [[`thrust::device_vector`, `thrust::host_vector`] [`thrust_algebra`] [`thrust_operations`] [For running odeint on CUDA devices by using __thrust]]
+  [[Any RandomAccessRange] [`openmp_range_algebra`] [`default_operations`] [OpenMP-parallelised range algebra]]
+  [[`openmp_state`] [`openmp_algebra`] [`default_operations`] [OpenMP-parallelised algebra for split data]]
+  [[`boost::array` or anything which allocates the elements in a C-like manner] [`vector_space_algebra`] [`mkl_operations`] [Using the __intel_mkl in odeint for maximum performance. Currently, only the RK4 stepper is supported.]]
+]
+
+[endsect]
+
+[section Example expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Vector operation] [`algebra.for_each3( y , x1 , x2 , Operations::scale_sum2< Value1 , Value2 >( a1 , a2 ) )`] [void] [Calculates ['*y* = a1 *x1* + a2 *x2*]]]
+]
+
+[endsect]
+
+[endsect]
\ No newline at end of file
diff --git a/doc/concepts/state_wrapper.qbk b/doc/concepts/state_wrapper.qbk
new file mode 100644
index 0000000..768671d
--- /dev/null
+++ b/doc/concepts/state_wrapper.qbk
@@ -0,0 +1,39 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright 2011 Mario Mulansky
+  Copyright 2012 Karsten Ahnert
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+[section State Wrapper]
+
+[heading Description]
+
+The `State Wrapper` concept describes the way odeint creates temporary state objects to store intermediate results within the stepper's `do_step` methods.
+
+[heading Notation]
+
+[variablelist
+  [[`State`] [A type that is the `state_type` of the ODE]]
+  [[`WrappedState`] [A type that is a model of State Wrapper for the state type `State`.]] 
+  [[`x`] [Object of type `State`]]
+  [[`w`] [Object of type `WrappedState`]]
+]
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Get resizeability] [`is_resizeable< State >`] [`boost::false_type` or `boost::true_type`] [Returns `boost::true_type` if the `State` is resizeable, `boost::false_type` otherwise.]]
+  [[Create `WrappedState` type] [`state_wrapper< State >`] [`WrappedState`] [Creates the type for a `WrappedState` for the state type `State`]]
+  [[Constructor] [`WrappedState()`] [`WrappedState`] [Constructs a state wrapper with an empty state]]
+  [[Copy Constructor] [`WrappedState( w )`] [`WrappedState`] [Constructs a state wrapper with a state of the same size as the state in `w`]]
+  [[Get state] [`w.m_v`] [`State`] [Returns the `State` object of this state wrapper.]]
+]
+
+[endsect]
\ No newline at end of file
diff --git a/doc/concepts/stepper.qbk b/doc/concepts/stepper.qbk
new file mode 100644
index 0000000..a23c027
--- /dev/null
+++ b/doc/concepts/stepper.qbk
@@ -0,0 +1,93 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright 2011-2012 Karsten Ahnert
+  Copyright 2011 Mario Mulansky
+  Copyright 2012 Sylwester Arabas
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+[section Stepper]
+
+This concepts specifies the interface a simple stepper has to fulfill to be used within the __integrate_functions.
+
+[heading Description]
+
+The basic stepper concept.
+A basic stepper following this Stepper concept is able to perform a single step of the solution /x(t)/ of an ODE to obtain /x(t+dt)/ using a given step size /dt/.
+Basic steppers can be Runge-Kutta steppers, symplectic steppers as well as implicit steppers.
+Depending on the actual stepper, the ODE is defined as __system, __symplectic_system, __simple_symplectic_system or __implicit_system.
+Note that all error steppers are also basic steppers.
+
+[heading Refinement of]
+
+* DefaultConstructable
+* CopyConstructable
+
+
+[heading Associated types]
+
+*  '''<para>'''[*state_type]'''</para>'''
+'''<para>'''`Stepper::state_type`'''</para>'''
+'''<para>'''The type characterizing the state of the ODE, hence ['x].'''</para>'''
+
+*  '''<para>'''[*deriv_type]'''</para>'''
+'''<para>'''`Stepper::deriv_type`'''</para>'''
+'''<para>'''The type characterizing the derivative of the ODE, hence ['d x/dt].'''</para>'''
+
+*  '''<para>'''[*time_type]'''</para>'''
+'''<para>'''`Stepper::time_type`'''</para>'''
+'''<para>'''The type characterizing the dependent variable of the ODE, hence the time ['t].'''</para>'''
+
+*  '''<para>'''[*value_type]'''</para>'''
+'''<para>'''`Stepper::value_type`'''</para>'''
+'''<para>'''The numerical data type which is used within the stepper, something like `float`, `double`, `complex&lt; double &gt;`.'''</para>'''
+
+*  '''<para>'''[*order_type]'''</para>'''
+'''<para>'''`Stepper::order_type`'''</para>'''
+'''<para>'''The type characterizing the order of the ODE, typically `unsigned short`.'''</para>'''
+
+*  '''<para>'''[*stepper_category]'''</para>'''
+'''<para>'''`Stepper::stepper_category`'''</para>'''
+'''<para>'''A tag type characterizing the category of the stepper. This type must be convertible to `stepper_tag`.'''</para>'''
+
+
+[heading Notation]
+
+[variablelist 
+  [[`Stepper`] [A type that is a model of Stepper]]
+  [[`State`] [A type representing the state /x/ of the ODE]]
+  [[`Time`] [A type representing the time /t/ of the ODE]]
+  [[`stepper`] [An object of type `Stepper`]]
+  [[`x`] [Object of type `State`]]
+  [[`t`, `dt`] [Objects of type `Time`]]
+  [[`sys`] [An object defining the ODE. Depending on the Stepper this might be a model of __system, __symplectic_system, __simple_symplectic_system or __implicit_system ]]
+]
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Get the order] [`stepper.order()`] [`order_type`] [Returns the order of the stepper.]]
+  [[Do step] [`stepper.do_step( sys , x , t , dt )`] [`void`] [Performs one step of step size `dt`.  The newly obtained state is written in place in `x`.] ]
+
+  [/ [Do step with reference] [`stepper.do_step( boost::ref(sys) , x , t , dt )`] [`void`] [Performs one step of step size `dt`.  The newly obtained state is written in place in `x`.] ]
+
+  [/ [Do step out-of-place] [`stepper.do_step( sys , in , t , out , dt )`] [`void`] [Performs one step. The newly obtained state is written to `out`] ]
+]
+
+[heading Models]
+
+* `runge_kutta4`
+* `euler`
+* `runge_kutta_cash_karp54`
+* `runge_kutta_dopri5`
+* `runge_kutta_fehlberg78`
+* `modified_midpoint`
+* `rosenbrock4`
+
+[endsect]
diff --git a/doc/concepts/symplectic_system.qbk b/doc/concepts/symplectic_system.qbk
new file mode 100644
index 0000000..e078597
--- /dev/null
+++ b/doc/concepts/symplectic_system.qbk
@@ -0,0 +1,99 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright 2011 Mario Mulansky
+  Copyright 2011-2012 Karsten Ahnert
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+[section Symplectic System]
+
+[heading Description]
+
+This concept describes how to define a symplectic system written with generalized coordinate `q` and generalized momentum `p`:
+
+[' q'(t) = f(p) ]
+
+[' p'(t) = g(q) ]
+
+Such a situation is typically found for Hamiltonian systems with a separable Hamiltonian:
+
+[' H(p,q) = H[sub kin](p) + V(q) ]
+
+which gives the equations of motion:
+
+[' q'(t) = dH[sub kin] / dp = f(p) ]
+
+[' p'(t) = dV / dq = g(q) ]
+
+
+The algorithmic implementation of this situation is described by a pair of callable objects for /f/ and /g/ with a specific parameter signature.
+Such a system should be implemented as a std::pair of functions or a functors.
+Symplectic systems are used in symplectic steppers like `symplectic_rkn_sb3a_mclachlan`.
+
+[heading Notation]
+
+[variablelist
+  [[`System`] [A type that is a model of SymplecticSystem]]
+  [[`Coor`] [The type of the coordinate  ['q]]]
+  [[`Momentum`] [The type of the momentum ['p]]]
+  [[`CoorDeriv`] [The type of the derivative of coordinate ['q']]]
+  [[`MomentumDeriv`] [The type of the derivative of momentum ['p']]]
+  [[`sys`] [An object of the type `System`]]
+  [[`q`] [Object of type Coor]]
+  [[`p`] [Object of type Momentum]]
+  [[`dqdt`] [Object of type CoorDeriv]]
+  [[`dpdt`] [Object of type MomentumDeriv]]
+]
+
+[heading Valid expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Check for pair] [`boost::is_pair< System >::type`] [`boost::mpl::true_`] [Check if System is a pair]]
+  [[Calculate ['dq/dt = f(p)]] [`sys.first( p , dqdt )`] [`void`] [Calculates ['f(p)], the result is stored into `dqdt`] ]
+  [[Calculate ['dp/dt = g(q)]] [`sys.second( q , dpdt )`] [`void`] [Calculates ['g(q)], the result is stored into `dpdt`] ]
+]
+
+[endsect]
+
+
+[section Simple Symplectic System]
+
+[heading Description]
+
+In most Hamiltonian systems the kinetic term is a quadratic term in the momentum ['H[sub kin] = p^2 / 2m] and in many cases it is possible to rescale coordinates and set /m=1/ which leads to a trivial equation of motion:
+
+[' q'(t) = f(p) = p. ]
+
+while for /p'/ we still have the general form
+
+[' p'(t) = g(q) ]
+
+As this case is very frequent we introduced a concept where only the nontrivial equation for /p'/ has to be provided to the symplectic stepper.
+We call this concept ['SimpleSymplecticSystem]
+
+[heading Notation]
+
+[variablelist
+  [[System] [A type that is a model of SimpleSymplecticSystem]]
+  [[Coor] [The type of the coordinate  ['q]]]
+  [[MomentumDeriv] [The type of the derivative of momentum ['p']]]
+  [[sys] [An object that models System]]
+  [[q] [Object of type Coor]]
+  [[dpdt] [Object of type MomentumDeriv]]
+]
+
+[heading Valid Expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Check for pair] [`boost::is_pair< System >::type`] [`boost::mpl::false_`] [Check if System is a pair, should be evaluated to false in this case.]]
+  [[Calculate ['dp/dt = g(q)]] [`sys( q , dpdt )`] [`void`] [Calculates ['g(q)], the result is stored into `dpdt`] ]
+]
+
+[endsect]
\ No newline at end of file
diff --git a/doc/concepts/system.qbk b/doc/concepts/system.qbk
new file mode 100644
index 0000000..d04c998
--- /dev/null
+++ b/doc/concepts/system.qbk
@@ -0,0 +1,43 @@
+[/============================================================================
+  Boost.odeint
+
+  Copyright 2011 Mario Mulansky
+  Copyright 2011-2012 Karsten Ahnert
+
+  Use, modification and distribution is subject to 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)
+=============================================================================/]
+
+
+[section System]
+
+[heading Description]
+
+The System concept models the algorithmic implementation of the rhs. of the ODE ['x' = f(x,t)].
+The only requirement for this concept is that it should be callable with a specific parameter syntax (see below).
+A System is typically implemented as a function or a functor.
+Systems fulfilling this concept are required by all Runge-Kutta steppers as well as the Bulirsch-Stoer steppers.
+However, symplectic and implicit steppers work with other system concepts, see __symplectic_system and __implicit_system.
+
+[heading Notation]
+
+[variablelist
+  [[`System`] [A type that is a model of System]] 
+  [[`State`] [A type representing the state /x/ of the ODE]]
+  [[`Deriv`] [A type representing the derivative /x'/ of the ODE]]
+  [[`Time`] [A type representing the time]]
+  [[`sys`] [An object of type `System`]]
+  [[`x`] [Object of type `State`]]
+  [[`dxdt`] [Object of type `Deriv`]]
+  [[`t`] [Object of type `Time`]]
+]
+
+[heading Valid expressions]
+
+[table
+  [[Name] [Expression] [Type] [Semantics]]
+  [[Calculate ['dx/dt := f(x,t)]] [`sys( x , dxdt , t )`] [`void`] [Calculates f(x,t), the result is stored into dxdt] ]
+]
+
+[endsect]
\ No newline at end of file