diff --git a/doc/sphinx/reactors.rst b/doc/sphinx/reactors.rst
index 103a412bf..07dff227f 100644
--- a/doc/sphinx/reactors.rst
+++ b/doc/sphinx/reactors.rst
@@ -1,74 +1,64 @@
 .. default-role:: math
 
-****************
-Reactor Networks
-****************
+.. py:currentmodule:: cantera
 
-Cantera's Reactor Network module is designed to simulate networks of
-interconnected reactors. The contents of each reactor in the network are
-assumed to be homogeneous, a model variously referred to as the Continuously
-Stirred Tank Reactor (CSTR), Well-Stirred Reactor (WSR), or Perfectly Stirred
-Reactor (PSR) model. Cantera solves the time-dependent governing equations
-that describe the evolution of the chemical and thermodynamic state of the
-reactors.
+*****************************
+Reactors and Reactor Networks
+*****************************
 
-The contents of each reactor can undergo chemical reactions according to a
-specified kinetic mechanism, and surface reactions may occur on the reactor
-walls. Each reactor in a network may be connected so that the contents of one
-reactor flow into another. Reactors may be also be in contact with one another
-or the environment via walls which move or conduct heat.
+A Cantera Reactor represents the simplest form of a chemically reacting system.
+It corresponds to an extensive thermodynamic control volume `V`, in which all
+state variables are homogeneously distributed. The system is generally unsteady,
+i.e. all states are functions of time. In particular, transient state changes
+due to chemical reactions are possible. However, thermodynamic (but not
+chemical) equilibrium is assumed to be present throughout the reactor at all
+instants of time.
 
-The purpose of this document is to describe the governing equations of reactor
-models as implemented in Cantera.
+Reactors can interact with the surrounding environment in multiple ways:
 
-Wall Interactions
-=================
+- Expansion/compression work: By moving the walls of the reactor, its volume can
+  be changed and expansion or compression work can be done by or on the system,
+  i.e., the Reactor.
+- Heat transfer: An arbitrary heat transfer rate can be defined to cross the
+  boundaries of the reactor.
+- Mass transfer: The reactor can have multiple inlets and outlets. For the
+  inlets, arbitrary states can be defined. Through the outlets, fluid with the
+  current state of the reactor exits the reactor.
+- Surface interaction: One or multiple walls can influence the chemical
+  reactions in the reactor. This is not just restricted to catalytic reactions,
+  but mass transfer between the surface and the fluid can also be modeled.
 
-At each wall where there are surface reactions, there is a net generation (or
-destruction) of homogeneous phase species. The molar rate of production for
-each species `k` on wall `w` is `\dot{s}_{k,w}`. The total (mass) production
-rate for species `k` on all walls is:
+All of these interactions do not have to be constant, but can vary as a function
+of time or state. For example, heat transfer can be described as a function of
+the temperature difference between the reactor and the environment, or the wall
+movement can be modeled depending on the pressure difference. Typically,
+interactions of the reactor with the environment are defined on one or multiple
+*walls*, *inlets*, and *outlets*.
 
-.. math::
+In addition to single reactors, Cantera is also able to interconnect reactors 
+into a *Reactor Network*. Each reactor in a network may be connected so that 
+the contents of one reactor flow into another. Reactors may also be in contact 
+with one another or the environment via walls which move or conduct heat.
 
-    \dot{m}_{k,wall} = W_k \sum_w A_w \dot{s}_{k,w}
-
-where `W_k` is the molecular weight of species `k` and `A_w` is the area of
-each wall. The net mass flux from all walls is then:
-
-.. math::
-
-    \dot{m}_{wall} = \sum_k \dot{m}_{k,wall}
-
-The total rate of heat transfer through all walls is:
-
-.. math::
-
-    \dot{Q} = \sum_w f_w \dot{Q}_w
-
-General Reactor
-===============
+Governing Equations for Single Reactors
+=======================================
 
 The state variables for Cantera's general reactor model are
 
-    - `m`, the mass of the reactor's contents
-    - `V`, the reactor volume
-    - `U`, the total internal energy of the reactors contents
-    - `Y_k`, the mass fractions for each species
+- `m`, the mass of the reactor's contents (in kg)
+- `V`, the reactor volume (in m\ :sup:`3`) (not a state variable for
+  *Constant Pressure Reactor* and *Ideal Gas Constant Pressure Reactor*)
+- A state variable describing the energy of the system, depending on the
+  configuration (see `Energy Conservation`_ for further explanation):
 
-Reactor Volume
---------------
+  - General *Reactor*: `U`, the total internal energy of the reactors
+    contents (in J)
+  - *Constant Pressure Reactor*: `H`, the total enthalpy of the reactors
+    contents (in J)
+  - *Ideal Gas Reactor* and *Ideal Gas Constant Pressure Reactor*: `T`, the
+    temperature (in K)
 
-The reactor volume changes as a function of time due to the motion of one or
-more walls:
-
-.. math::
-
-    \frac{dV}{dt} = \sum_w f_w A_w v_w(t)
-
-where `f_w = \pm 1` indicates the facing of the wall, `A_w` is the surface
-area of the wall, and `v_w(t)` is the velocity of the wall as a function of
-time.
+- `Y_k`, the mass fractions for each species (dimensionless)
 
 Mass Conservation
 -----------------
@@ -86,8 +76,8 @@ Species Conservation
 --------------------
 
 The rate at which species `k` is generated through homogeneous phase reactions
-is `V \dot{\omega}_k W_k`, and the total rate at which species `k` is
-generated is:
+is `V \dot{\omega}_k W_k`, and the total rate at which species `k` is generated
+is:
 
 .. math::
 
@@ -106,12 +96,40 @@ for `dm/dt`, the equation for each homogeneous phase species is:
 
 .. math::
 
-    m \frac{dY}{dt} = \sum_{in} \dot{m}_{in} (Y_{k,in} - Y_k)+
+    m \frac{dY_k}{dt} = \sum_{in} \dot{m}_{in} (Y_{k,in} - Y_k)+
                       \dot{m}_{k,gen} - Y_k \dot{m}_{wall}
 
+
+Reactor Volume
+--------------
+
+The reactor volume changes as a function of time due to the motion of one or
+more walls:
+
+.. math::
+
+    \frac{dV}{dt} = \sum_w f_w A_w v_w(t)
+
+where `f_w = \pm 1` indicates the facing of the wall, `A_w` is the surface
+area of the wall, and `v_w(t)` is the velocity of the wall as a function of
+time.
+
+For *Constant Pressure Reactor* and *Ideal Gas Constant Pressure Reactor*, the
+volume is not a state variable, but instead takes on whatever value is
+consistent with holding the pressure constant.
+
 Energy Conservation
 -------------------
 
+The solution of the energy equation can be enabled or disabled by changing the 
+``energy_enabled`` flag. It is enabled by default.
+
+The implemented formulation of the energy equation depends on which reactor 
+model is used.
+
+Standard Reactor
+****************
+
 The equation for the total internal energy is found by writing the first law
 for an open system:
 
@@ -120,13 +138,36 @@ for an open system:
     \frac{dU}{dt} = - p \frac{dV}{dt} - \dot{Q} +
                     \sum_{in} \dot{m}_{in} h_{in} - h \sum_{out} \dot{m}_{out}
 
-Ideal Gas Reactor
-=================
 
-The Ideal Gas Reactor model is similar to the General Reactor model, with the
-reactor temperature `T` replacing the total internal energy `U` as a state
-variable. For an ideal gas, we can rewrite the total internal energy in terms
-of the mass fractions and temperature:
+Constant Pressure Reactor
+*************************
+
+For this reactor model, the pressure is held constant. The volume is not a
+state variable, but instead takes on whatever value is consistent with holding
+the pressure constant. The total enthalpy replaces the total internal energy
+as a state variable. Using the definition of the total enthalpy:
+
+.. math::
+
+    H = U + pV
+
+    \frac{d H}{d t} = \frac{d U}{d t} + p \frac{dV}{dt} + V \frac{dp}{dt}
+
+Noting that `dp/dt = 0` and substituting into the energy equation yields:
+
+.. math::
+
+   \frac{dH}{dt} = - \dot{Q} + \sum_{in} \dot{m}_{in} h_{in}
+                   - h \sum_{out} \dot{m}_{out}
+
+
+Ideal Gas Reactor
+*****************
+
+In case of the Ideal Gas Reactor Model, the reactor temperature `T` is used 
+instead of the total internal energy `U` as a state variable. For an ideal gas, 
+we can rewrite the total internal energy in terms of the mass fractions and 
+temperature:
 
 .. math::
 
@@ -150,34 +191,11 @@ significantly reduces the cost of evaluating the system Jacobian, since the
 derivatives of the species equations are taken at constant temperature instead
 of constant internal energy.
 
-Constant Pressure Reactor
-=========================
-
-For this reactor model, the pressure is held constant. The volume is not a
-state variable, but instead takes on whatever value is consistent with holding
-the pressure constant. The total enthalpy replaces the total internal energy
-as a state variable. Using the definition of the total enthalpy:
-
-.. math::
-
-    H = U + pV
-
-    \frac{dH}{dt} = p \frac{dV}{dt} + V \frac{dp}{dt}
-
-Noting that `dp/dt = 0` and substituting into the energy equation yields:
-
-.. math::
-
-   \frac{dH}{dt} = - \dot{Q} + \sum_{in} \dot{m}_{in} h_{in}
-                   - h \sum_{out} \dot{m}_{out}
-
-The species and continuity equations are the same as for the general reactor
-model.
 
 Ideal Gas Constant Pressure Reactor
-===================================
+***********************************
 
-As for the Ideal Gas Reactor, we replace the total enthalpy as a state
+As for the Ideal Gas Reactors, we replace the total enthalpy as a state
 variable with the temperature by writing the total enthalpy in terms of the
 mass fractions and temperature:
 
@@ -195,3 +213,378 @@ temperature:
 
     m c_p \frac{dT}{dt} = - \dot{Q} - \sum_k h_k \dot{m}_{k,gen}
         + \sum_{in} \dot{m}_{in} \left(h_{in} - \sum_k h_k Y_{k,in} \right)
+
+
+Wall Interactions
+-----------------
+
+The total rate of heat transfer through all walls is:
+
+.. math::
+
+    \dot{Q} = \sum_w f_w \dot{Q}_w
+
+where `f_w = \pm 1` indicates the facing of the wall.
+
+In case of surface reactions, there is a net generation (or
+destruction) of homogeneous phase species at the wall. The molar rate of 
+production for each species `k` on wall `w` is `\dot{s}_{k,w}` (in kmol/s/m\ 
+:sup:`2`). The total (mass) production rate for species `k` on all walls is:
+
+.. math::
+
+    \dot{m}_{k,wall} = W_k \sum_w A_w \dot{s}_{k,w}
+
+where `W_k` is the molecular weight of species `k` and `A_w` is the area of
+each wall. The net mass flux from all walls is then:
+
+.. math::
+
+    \dot{m}_{wall} = \sum_k \dot{m}_{k,wall}
+
+
+Reactor Networks and Devices
+============================
+
+While reactors by themselves just define the above governing equations of the 
+reactor, the time integration is performed in reactor networks. A reactor 
+network is therefore necessary even if only a single reactor is considered.
+
+The advantage of reactor networks obviously is that multiple reactors can be 
+interconnected. Not only mass flow from one reactor into another can be 
+realized, but also heat can be transferred, or the wall between reactors can 
+move. To set up a network, the following components can be defined in addition 
+to the reactors previously mentioned:
+
+- **Reservoir**: A reservoir can be thought of as an infinitely large volume, in
+  which all states are predefined and never change from their initial values.
+  Typically, it represents a vessel to define temperature and composition of a
+  stream of mass flowing into a reactor, or the ambient fluid surrounding the
+  reactor network. Besides, the fluid flow finally finally exiting a reactor
+  network has to flow into a reservoir. In the latter case, the state of the
+  reservoir (except pressure) is irrelevant.
+
+- **Wall**: A wall separates two reactors, or a reactor and a reservoir. A wall
+  has a finite area, may conduct or radiate heat between the two reactors on
+  either side, and may move like a piston.
+
+  Walls are stateless objects in Cantera, meaning that no differential equation
+  is integrated to determine any wall property. Since it is the wall (piston)
+  velocity that enters the energy equation, this means that it is the velocity,
+  not the acceleration or displacement, that is specified. The wall velocity is
+  computed from
+
+  .. math:: v = K(P_{\rm left} - P_{\rm right}) + v_0(t),
+
+  where :math:`K` is a non-negative constant, and :math:`v_0(t)` is a specified
+  function of time. The velocity is positive if the wall is moving to the right.
+
+  The heat flux through the wall is computed from
+
+  .. math:: q = U(T_{\rm left} - T_{\rm right}) + \epsilon\sigma (T_{\rm left}^4 
+                - T_{\rm right}^4) + q_0(t),
+
+  where :math:`U` is the overall heat transfer coefficient for
+  conduction/convection, and :math:`\epsilon` is the emissivity. The function
+  :math:`q_0(t)` is a specified function of time. The heat flux is positive when 
+  heat flows from the reactor on the left to the reactor on the right.
+
+  A heterogeneous reaction mechanism may be specified for one or both of the
+  wall surfaces. The mechanism object (typically an instance of class Interface)
+  must be constructed so that it is properly linked to the object representing
+  the fluid in the reactor the surface in question faces. The surface
+  temperature on each side is taken to be equal to the temperature of the
+  reactor it faces.
+
+  Source: `Python <cython/zerodim.html#wall>`_ | :ct:`C++ <Wall>`
+
+- **Valve**: A valve is a flow devices with mass flow rate that is a function of
+  the pressure drop across it. The default behavior is linear:
+
+  .. math:: \dot m = K_v (P_1 - P_2)
+
+  if :math:`P_1 > P_2.` Otherwise, :math:`\dot m = 0`. However, an arbitrary
+  function can also be specified, such that
+
+  .. math:: \dot m = F(P_1 - P_2)
+
+  if :math:`P_1 > P_2`, or :math:`\dot m = 0` otherwise. It is never possible
+  for the flow to reverse and go from the downstream to the upstream
+  reactor/reservoir through a line containing a Valve object.
+
+  Valve objects are often used between an upstream reactor and a downstream
+  reactor or reservoir to maintain them both at nearly the same pressure. By
+  setting the constant :math:`K_v` to a sufficiently large value, very small
+  pressure differences will result in flow between the reactors that counteracts
+  the pressure difference.
+
+- **Mass Flow Controller**: A mass flow controller maintains a specified mass
+  flow rate independent of upstream and downstream conditions. The equation used
+  to compute the mass flow rate is
+
+  .. math:: \dot m = \max(\dot m_0, 0.0)
+
+  where :math:`\dot m_0` is either a constant value or a function of time. Note 
+  that if :math:`\dot m_0 < 0`, the mass flow rate will be set to zero, since 
+  reversal of the flow direction is not allowed.
+
+  Unlike a real mass flow controller, a MassFlowController object will maintain
+  the flow even if the downstream pressure is greater than the upstream
+  pressure. This allows simple implementation of loops, in which exhaust gas
+  from a reactor is fed back into it through an inlet. But note that this
+  capability should be used with caution, since no account is taken of the work
+  required to do this.
+
+- **Pressure Controller**: A pressure controller is designed to be used in
+  conjunction with another 'master' flow controller, typically a
+  MassFlowController. The master flow controller is installed on the inlet of
+  the reactor, and the corresponding PressureController is installed on on
+  outlet of the reactor. The PressureController mass flow rate is equal to the
+  master mass flow rate, plus a small correction dependent on the pressure
+  difference:
+
+  .. math:: \dot m = \dot m_{\rm master} + K_v(P_1 - P_2).
+
+Time Integration
+----------------
+
+Cantera provides an ODE solver for solving the stiff equations of reacting 
+systems. If installed in combination with SUNDIALS, their optimized solver is 
+used. Starting off the current state of the system, it can be advanced in time 
+by two methods:
+
+- ``step()``: The step method computes the state of the system at the a priori
+  unspecified time `t_{\rm new}`. The time `t_{\rm new}` is internally computed
+  so that all states of the system only change within a (specifiable) band of
+  absolute and relative tolerances. Additionally, the time step must not be
+  larger than a predefined maximum time step `\Delta t_{\rm max}`. The new time
+  `t_{\rm new}` is returned by this function.
+
+- ``advance``\ `(t_{\rm new})`: This method computes the state of the system at
+  time `t_{\rm new}`. `t_{\rm new}` describes the absolute time from the initial
+  time of the system. By calling this method in a for loop for pre-defined
+  times, the state of the system is obtained for exactly the times specified.
+  Internally, several ``step()`` calls are typically performed to reach the
+  accurate state at time `t_{\rm new}`.
+
+The use of the ``advance`` method in a loop has the advantage that it produces 
+results corresponding to a predefined time series. These are associated with a 
+predefined memory consumption and well comparable between simulation runs with 
+different parameters. However, some detail (e.g. a fast ignition process) might 
+not be resolved in the output data due to the typically large time steps.
+
+The ``step`` method results in much more data points because of the small 
+timesteps needed. Additionally, the absolute time has to be kept tracked of 
+manually.
+
+Even though Cantera comes pre-defined with typical parameters for tolerances 
+and the maximum internal time step, the solution sometimes diverges. To solve 
+this problem, three parameters can be tuned: The absolute time stepping 
+tolerances, the relative time stepping tolerances, and the maximum time step. A 
+reduction of the latter value is particularly useful when dealing with abrupt 
+changes in the boundary conditions (e.g. opening/closing valves, see also 
+example :ref:`py-example-ic_engine.py`).
+
+General Usage in Cantera
+========================
+
+In Cantera, the following steps are typically necessary to investigate a 
+reactor network:
+
+1. Define ``Solution`` objects for the fluids to be flowing through your 
+reactor network.
+
+2. Define the reactor type(s) and reservoir(s) that describe your system. Chose 
+Ideal Gas (Constant Pressure) Reactor(s) if you only consider ideal gas phases.
+
+3. *Optional:* Set up the boundary conditions and flow devices between reactors 
+or reservoirs.
+
+4. Define a reactor network which contains all the reactors previously created.
+
+5. Advance the simulation in time, typically in a for- or while-loop. Note that 
+only the current state is stored in Cantera by default. If you want to observe 
+the transient states, you manually have to keep track of them.
+
+6. Analyze the data.
+
+Note that Cantera always solves a transient problem. If you are interested in 
+steady-state conditions, you can run your simulation for a long time until the 
+states are converged (see e.g. example :ref:`py-example-surf_pfr.py`,
+:ref:`py-example-combustor.py`).
+
+Cantera comes with a broad variety of well-commented example scrips for reactor 
+networks. Please refer to them for further information (:ref:`Python <sec-cython-examples>`, :ref:`Matlab <sec-matlab-examples>`).
+
+Common Reactor Types and their Implementation in Cantera
+========================================================
+
+Batch Reactor at Constant Volume or at Constant Pressure
+--------------------------------------------------------
+
+If you are interested in how a homogeneous chemical composition changes in time 
+when it is left to its own, a simple batch reactor can be used. Two versions 
+are commonly considered: A rigid vessel with fixed volume but variable 
+pressure, or a system idealized at constant pressure but varying volume.
+
+In Cantera, such a simulation can be performed very easily. The initial state 
+of the solution can be specified by composition and a set of thermodynamic 
+parameters (like temperature and pressure) as a standard Cantera solution 
+object. Upon its base, a general (Ideal Gas) Reactor or an (Ideal Gas) Constant 
+Pressure Reactor can be created, depending on if a constant volume or constant 
+pressure batch reactor should be considered, respectively. The behavior of the 
+solution in time can be simulated as a very simple Reactor Network containing 
+only the formerly created reactor.
+
+An example for such a Batch Reactor is :ref:`py-example-reactor1.py`.
+
+Continuously Stirred Tank Reactor
+---------------------------------
+
+A Continuously Stirred Tank Reactor (CSTR), also often referred to as 
+Well-Stirred Reactor (WSR), Perfectly Stirred Reactor (PSR), or Longwell 
+Reactor, is essentially a single Cantera reactor with an inlet, an outlet, and 
+constant volume. Therefore, the `Governing Equations for Single Reactors`_ 
+defined above apply accordingly.
+
+Steady state solutions to CSTRs are often of interest. In this case, the mass 
+flow rate `\dot{m}` is constant and equal at inlet and outlet. The mass 
+contained in the confinement `m` divided by `\dot{m}` defines the mean 
+residence time of the fluid in the confinement.
+
+At steady state, the time derivatives in the governing equations become zero, 
+and the system of ordinary differential equations can be reduced to a set of 
+coupled nonlinear algebraic equations. A Newton solver could be used to solve 
+this system of equations. However, a sophisticated implementation might be 
+required to account for the strong nonlinearities and the presence of multiple 
+solutions.
+
+Cantera does not have such a Newton solver implemented. Instead, steady CSTRs 
+are simulated by considering a time-dependent constant volume reactor with 
+specified in- and outflow conditions. Starting off at an initial solution, the 
+reactor network containing this reactor is advanced in time until the state of 
+the solution is converged. An example for this procedure is
+:ref:`py-example-combustor.py`.
+
+A problem can be the ignition of a CSTR: If the reactants are not reactive 
+enough, the simulation can result in the trivial solution that inflow and 
+outflow states are identical. To solve this problem, the reactor can be 
+initialized with a high temperature and/or radical concentration. A good 
+approach is to use the equilibrium composition of the reactants (which can be 
+computed using Cantera's ``equilibrate`` function) as an initial guess.
+
+
+Plug-Flow Reactor
+-----------------
+
+A Plug-Flow Reactor (PFR) represents a steady-state channel with a 
+cross-sectional area `A`. Typically an ideal gas flows through it at a constant 
+mass flow rate `\dot{m}`. Perpendicular to the flow direction, the gas is 
+considered to be completely homogeneous. In the axial direction `z`, the states 
+of the gas is allowed to change. However, all diffusion processes are neglected.
+
+Plug-Flow Reactors are often used to simulate ignition delay times, emission 
+formation, and catalytic processes.
+
+The governing equations of Plug-Flow Reactors are [KCG2003]_:
+
+- Mass conservation:
+
+  .. math:: \frac{d(\rho u A)}{dz} =  P' \sum_k \dot{s}_k W_k
+
+  where `u` is the axial velocity in (m/s) and `P'` is the chemically active 
+  channel perimeter in (m) (chemically active perimeter per unit length).
+
+- Continuity equation of species `k`:
+
+ .. math:: \rho u \frac{d Y_k}{dz} + Y_k P' \sum_k \dot{s}_k W_k =
+      \dot{\omega}_k W_k + P' \dot{s}_k W_k
+
+- Energy conservation:
+
+  .. math:: \rho u A c_p \frac{d T}{d z} =
+     - A \sum_k h_k \dot{\omega}_k W_k
+     - P' \sum_k h_k \dot{s}_k W_k
+     + U P (T_w - T)
+
+  where `U` is the heat transfer coefficient in (W/m/K), `P` is the perimeter of 
+  the duct in (m), and `T_w` is the wall temperature in (K). Kinetic and 
+  potential energies are neglected.
+
+- Momentum conservation in the axial direction:
+
+  .. math:: \rho u A \frac{d u}{d z} + u P' \sum_k \dot{s}_k W_k =
+     - \frac{d (p A)}{dz} - \tau_w P
+
+  where `\tau_w` is the wall friction coefficient (which might be computed from 
+  Reynolds number based correlations).
+
+Even though this problem extends geometrically in one direction, it can be 
+modeled via zero-dimensional reactors: Due to the neglecting of diffusion, 
+downstream parts of the reactor have no influence on upstream parts. Therefore, 
+PFRs can be modeled by marching from the beginning to the end of the reactor.
+
+Cantera does not (yet) provide dedicated class to solve the PFR equations (The 
+``FlowReactor`` class is currently under development). However, there are two 
+ways to simulate a PFR with the reactor elements previously presented. Both 
+rely on the assumption that pressure is approximately constant throughout the 
+Plug-Flow Reactor and that there is no friction. The momentum conservation 
+equation is thus neglected.
+
+
+PFR Modeling by Considering a Lagrangian Reactor
+************************************************
+
+A Plug-Flow Reactor can also be described from a Lagrangian point of view: An 
+unsteady fluid particle is considered which travels along the axial streamline 
+through the PFR. Since there is no information traveling upstream, the state 
+change of the fluid particle can be computed by a forward (upwind) integration 
+in time. Using the continuity equation, the speed of the particle can be 
+derived. By integrating the velocity in time, the temporal information can be 
+translated into the spatial resolution of the PFR.
+
+An example for this procedure can be found in :ref:`py-example-pfr.py`.
+
+
+PFR Modeling as a Series of CSTRs
+*********************************
+
+The Plug-Flow Reactor is spatially discretized into a large number of axially 
+distributed volumes. These volumes are modeled to be steady-state CSTRs.
+
+The only reason to use this approach as opposed to the Lagrangian one is if you 
+need to include surface reactions, because the system of equations ends up 
+being a DAE system instead of an ODE system.
+
+In Cantera, it is sufficient to consider a single reactor and march it forward 
+in time, because there is no information traveling upstream. The mass flow rate 
+`\dot{m}` through the PFR enters the reactor from an upstream reservoir. For 
+the first reactor, the reservoir conditions are the inflow boundary conditions 
+of the PFR. By performing a time integration as described in `Continuously 
+Stirred Tank Reactor`_ until the state of the reactor is converged, the 
+steady-state CSTR solution is computed. The state of the CSTR is the inlet 
+boundary condition for the next CSTR downstream.
+
+An example for this procedure can be found in :ref:`py-example-pfr.py` and
+:ref:`py-example-surf_pfr.py`.
+
+
+Advanced Concepts
+=================
+
+In some cases, Cantera's solver is insufficient to describe a certain 
+configuration. In this situation, Cantera can still be used to provide chemical 
+and thermodynamic computations, but external ODE solvers can be applied. See 
+example :ref:`py-example-custom.py`.
+
+
+Literature
+==========
+
+For further reading, the following books are recommended:
+
+.. [KCG2003] Kee, Coltrin, Glarborg: *Chemically Reacting Flow*. 
+             Wiley-Interscience, 2003
+
+.. [Tur2000] Turns: *An Introduction to Combustion: Concepts and Applications*, 
+             McGraw Hill, 2000
diff --git a/interfaces/cython/cantera/examples/reactors/ic_engine.py b/interfaces/cython/cantera/examples/reactors/ic_engine.py
new file mode 100644
index 000000000..6ad369549
--- /dev/null
+++ b/interfaces/cython/cantera/examples/reactors/ic_engine.py
@@ -0,0 +1,263 @@
+# -*- coding: utf-8 -*-
+"""
+Simulation of a (gaseous) Diesel-type internal combustion engine.
+
+The use of pure propane as fuel requires an unrealistically high compression
+ratio.
+
+"""
+
+import cantera as ct
+import numpy as np
+
+#######################################################################
+# Input Parameters
+#######################################################################
+
+f = 3000. / 60.  # engine speed [1/s] (3000 rpm)
+V_H = .5e-3  # displaced volume [m**3]
+epsilon = 50.  # compression ratio [-]
+d_piston = 0.083  # piston diameter [m]
+
+# turbocharger temperature, pressure, and composition
+T_inlet = 300.  # K
+p_inlet = 1.3e5  # Pa
+comp_inlet = 'O2:1, N2:3.76'
+
+# outlet pressure
+p_outlet = 1.2e5  # Pa
+
+# fuel properties (gaseous!)
+T_injector = 300.  # K
+p_injector = 1600e5  # Pa
+comp_injector = 'C3H8:1'
+
+# ambient properties
+T_ambient = 300.  # K
+p_ambient = 1e5  # Pa
+comp_ambient = 'O2:1, N2:3.76'
+
+# Reaction mechanism name
+reaction_mechanism = 'gri30.xml'
+
+# Inlet valve friction coefficient, open and close timings
+inlet_valve_coeff = 1.e-6
+inlet_open = -18. / 180. * np.pi
+inlet_close = 198. / 180. * np.pi
+
+# Outlet valve friction coefficient, open and close timings
+outlet_valve_coeff = 1.e-6
+outlet_open = 522. / 180 * np.pi
+outlet_close = 18. / 180. * np.pi
+
+# Fuel mass, injector open and close timings
+injector_open = 350. / 180. * np.pi
+injector_close = 365. / 180. * np.pi
+injector_mass = 3.2e-5  # kg
+injector_t_open = (injector_close - injector_open) / 2. / np.pi / f
+
+# Simulation time and resolution
+sim_n_revolutions = 8.
+sim_n_timesteps = 100000.
+
+###################################################################
+
+# load reaction mechanism
+gas = ct.Solution(reaction_mechanism)
+
+# define initial state
+gas.TPX = T_inlet, p_inlet, comp_inlet
+r = ct.IdealGasReactor(gas)
+# define inlet state
+gas.TPX = T_inlet, p_inlet, comp_inlet
+inlet = ct.Reservoir(gas)
+# define injector state (gaseous!)
+gas.TPX = T_injector, p_injector, comp_injector
+injector = ct.Reservoir(gas)
+# define outlet pressure (temperature and composition don't matter)
+gas.TPX = T_ambient, p_outlet, comp_ambient
+outlet = ct.Reservoir(gas)
+# define ambient pressure (temperature and composition don't matter)
+gas.TPX = T_ambient, p_ambient, comp_ambient
+ambient_air = ct.Reservoir(gas)
+
+# set up connecting devices
+inlet_valve = ct.Valve(inlet, r)
+injector_mfc = ct.MassFlowController(injector, r)
+outlet_valve = ct.Valve(r, outlet)
+piston = ct.Wall(ambient_air, r)
+
+# convert time to crank angle
+def crank_angle(t):
+    return np.remainder(2 * np.pi * f * t, 4 * np.pi)
+
+# set up IC engine parameters
+V_oT = V_H / (epsilon - 1.)
+A_piston = .25 * np.pi * d_piston ** 2
+stroke = V_H / A_piston
+r.volume = V_oT
+piston.area = A_piston
+def piston_speed(t):
+    return - stroke / 2 * 2 * np.pi * f * np.sin(crank_angle(t))
+piston.set_velocity(piston_speed)
+
+# create a reactor network containing the cylinder
+sim = ct.ReactorNet([r])
+
+# set up output data arrays
+t_sim = sim_n_revolutions / f
+t = (np.arange(sim_n_timesteps) + 1) / sim_n_timesteps * t_sim
+p = np.zeros_like(t)
+V = np.zeros_like(t)
+T = np.zeros_like(t)
+s = np.zeros_like(t)
+m = np.zeros_like(t)
+test = np.zeros_like(t)
+mdot_in = np.zeros_like(t)
+mdot_out = np.zeros_like(t)
+MW = np.zeros_like(t)
+d_W_v_d_t = np.zeros_like(t)
+heat_release_rate = np.zeros_like(t)
+species_X = np.zeros((t.size, gas.n_species))
+
+# set parameters for the automatic time step refinement
+n_last_refinement = -np.inf  # for initialization only
+n_wait_coarsening = 10
+
+# do simulation
+for n1, t_i in enumerate(t):
+    # define opening and closing of valves and injector
+    if (np.mod(crank_angle(t_i) - inlet_open, 4 * np.pi) <
+            np.mod(inlet_close - inlet_open, 4 * np.pi)):
+        inlet_valve.set_valve_coeff(inlet_valve_coeff)
+        test[n1] = 1
+    else:
+        inlet_valve.set_valve_coeff(0)
+    if (np.mod(crank_angle(t_i) - outlet_open, 4 * np.pi) <
+            np.mod(outlet_close - outlet_open, 4 * np.pi)):
+        outlet_valve.set_valve_coeff(outlet_valve_coeff)
+    else:
+        outlet_valve.set_valve_coeff(0)
+    if (np.mod(crank_angle(t_i) - injector_open, 4 * np.pi) <
+            np.mod(injector_close - injector_open, 4 * np.pi)):
+        injector_mfc.set_mass_flow_rate(injector_mass / injector_t_open)
+    else:
+        injector_mfc.set_mass_flow_rate(0)
+
+    # perform time integration, refine time step if necessary
+    for n2 in range(4):
+        if n2 is 4:
+            raise 'Error: Refinement limit reached'
+        try:
+            sim.advance(t_i)
+        except Exception:
+            sim.set_max_time_step(1e-6 * 10. ** -n2)
+            n_last_refinement = n1
+    # coarsen time step if too long ago
+    if n1 - n_last_refinement is n_wait_coarsening:
+        sim.set_max_time_step(1e-5)
+
+    # write output data
+    p[n1] = r.thermo.P
+    V[n1] = r.volume
+    T[n1] = r.T
+    s[n1] = r.thermo.s
+    m[n1] = r.mass
+    mdot_in[n1] = inlet_valve.mdot(0)
+    mdot_out[n1] = outlet_valve.mdot(0)
+    MW[n1] = r.thermo.mean_molecular_weight
+    species_X[n1] = r.thermo.X
+    d_W_v_d_t[n1] = - (r.thermo.P - ambient_air.thermo.P) * A_piston * \
+        piston_speed(t_i)
+    heat_release_rate[n1] = - r.volume * ct.gas_constant * r.T * \
+        np.sum(gas.standard_enthalpies_RT * r.thermo.net_production_rates, 0)
+
+
+#####################################################################
+# Plot Results in matplotlib
+#####################################################################
+
+import matplotlib.pyplot as plt
+
+# pressure and temperature
+plt.figure()
+plt.clf()
+plt.subplot(211)
+plt.plot(t, p / 1.e5)
+plt.ylabel('$p$ [bar]')
+plt.xlabel('$\phi$ [deg]')
+plt.xticks(plt.xticks()[0], [])
+plt.subplot(212)
+plt.plot(t, T)
+plt.ylabel('$T$ [K]')
+plt.xlabel('$\phi$ [deg]')
+plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
+           rotation=17)
+plt.show()
+plt.savefig('ic_engine_t_p_T.png')
+
+# p-V diagram
+plt.figure()
+plt.clf()
+plt.plot(V[t > 0.04] * 1000, p[t > 0.04] / 1.e5)
+plt.xlabel('$V$ [l]')
+plt.ylabel('$p$ [bar]')
+plt.show()
+plt.savefig('ic_engine_p_V.png')
+
+# T-S diagram
+plt.figure()
+plt.clf()
+plt.plot(m[t > 0.04] * s[t > 0.04], T[t > 0.04])
+plt.xlabel('$S$ [J/K]')
+plt.ylabel('$T$ [K]')
+plt.show()
+plt.savefig('ic_engine_T_S.png')
+
+# heat of reaction and expansion work
+plt.figure()
+plt.clf()
+plt.plot(t, heat_release_rate, label='$\dot{Q}$')
+plt.plot(t, d_W_v_d_t, label='$\dot{W}_v$')
+plt.ylim(-1e5, 1e6)
+plt.legend(loc=0)
+plt.ylabel('[W]')
+plt.xlabel('$\phi$ [deg]')
+plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
+           rotation=17)
+plt.show()
+plt.savefig('ic_engine_Q_W.png')
+
+# gas composition
+plt.figure()
+plt.clf()
+plt.plot(t, species_X[:, gas.species_index('O2')], label='O2')
+plt.plot(t, species_X[:, gas.species_index('CO2')], label='CO2')
+plt.plot(t, species_X[:, gas.species_index('CO')], label='CO')
+plt.plot(t, species_X[:, gas.species_index('C3H8')] * 10, label='C3H8 x10')
+plt.legend(loc=0)
+plt.ylabel('$X_i$ [-]')
+plt.xlabel('$\phi$ [deg]')
+plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
+           rotation=17)
+plt.show()
+plt.savefig('ic_engine_t_X.png')
+
+
+#####################################################################
+# Integral Results
+#####################################################################
+
+from scipy.integrate import trapz
+Q = trapz(heat_release_rate, t)
+W = trapz(d_W_v_d_t, t)
+eta = W / Q
+CO_emission = trapz(MW * mdot_out * species_X[:, gas.species_index('CO')], t) \
+    / trapz(MW * mdot_out, t)
+print('Heat release rate per cylinder (estimate):\t' +
+      format(Q / t_sim / 1000., ' 2.1f') + ' kW')
+print('Expansion power per cylinder (estimate):\t' +
+      format(W / t_sim / 1000., ' 2.1f') + ' kW')
+print('Efficiency (estimate):\t\t\t' + format(eta * 100., ' 2.1f') + ' %')
+print('CO emission (estimate):\t\t' + format(CO_emission * 1.e6, ' 2.1f') +
+      ' ppm')
diff --git a/interfaces/cython/cantera/examples/reactors/pfr.py b/interfaces/cython/cantera/examples/reactors/pfr.py
new file mode 100644
index 000000000..75d96fcdd
--- /dev/null
+++ b/interfaces/cython/cantera/examples/reactors/pfr.py
@@ -0,0 +1,174 @@
+# -*- coding: utf-8 -*-
+"""
+This example solves a plug-flow reactor problem of hydrogen-oxygen combustion.
+The PFR is computed by two approaches: The simulation of a Lagrangian fluid
+particle, and the simulation of a chain of reactors.
+"""
+
+import cantera as ct
+import numpy as np
+
+#######################################################################
+# Input Parameters
+#######################################################################
+
+T_0 = 1500.0  # inlet temperature [K]
+pressure = ct.one_atm  # constant pressure [Pa]
+composition_0 = 'H2:2, O2:1, AR:0.1'
+length = 1.5e-7  # *approximate* PFR length [m]
+u_0 = .006  # inflow velocity [m/s]
+area = 1.e-4  # cross-sectional area [m**2]
+
+# input file containing the reaction mechanism
+reaction_mechanism = 'h2o2.xml'
+
+# Resolution: The PFR will be simulated by 'n_steps' time steps or by a chain
+# of 'n_steps' stirred reactors.
+n_steps = 2000
+#####################################################################
+
+
+#####################################################################
+# Method 1: Lagrangian Particle Simulation
+#####################################################################
+# A Lagrangian particle is considered which travels through the PFR. Its
+# state change is computed by upwind time stepping. The PFR result is produced
+# by transforming the temporal resolution into spatial locations.
+# The spatial discretization is therefore not provided a priori but is instead
+# a result of the transformation.
+
+# import the gas model and set the initial conditions
+gas1 = ct.Solution(reaction_mechanism)
+gas1.TPX = T_0, pressure, composition_0
+mass_flow_rate1 = u_0 * gas1.density * area
+
+# create a new reactor
+r1 = ct.IdealGasConstPressureReactor(gas1)
+# create a reactor network for performing time integration
+sim1 = ct.ReactorNet([r1])
+
+# approximate a time step to achieve a similar resolution as in the next method
+t_total = length / u_0
+dt = t_total / n_steps
+# define time, space, and other information vectors
+t1 = (np.arange(n_steps) + 1) * dt
+z1 = np.zeros_like(t1)
+u1 = np.zeros_like(t1)
+T1 = np.zeros_like(t1)
+X_H2_1 = np.zeros_like(t1)
+for n1, t_i in enumerate(t1):
+    # perform time integration
+    sim1.advance(t_i)
+    # compute velocity and transform into space
+    u1[n1] = mass_flow_rate1 / area / r1.thermo.density
+    z1[n1] = z1[n1 - 1] + u1[n1] * dt
+    # write output data
+    T1[n1] = r1.T
+    X_H2_1[n1] = r1.thermo['H2'].X
+#####################################################################
+
+
+#####################################################################
+# Method 2: Chain of Reactors
+#####################################################################
+# The plug flow reactor is represented by a linear chain of zero-dimensional
+# reactors. The gas at the inlet to the first one has the specified inlet
+# composition, and for all others the inlet composition is fixed at the
+# composition of the reactor immediately upstream. Since in a PFR model there
+# is no diffusion, the upstream reactors are not affected by any downstream
+# reactors, and therefore the problem may be solved by simply marching from
+# the first to last reactor, integrating each one to steady state.
+# (This approach is anologous to the one presented in 'surf_pfr.py', which
+# additionally includes surface chemistry)
+
+
+# import the gas model and set the initial conditions
+gas2 = ct.Solution(reaction_mechanism)
+gas2.TPX = T_0, pressure, composition_0
+mass_flow_rate2 = u_0 * gas2.density * area
+dz = length / n_steps
+r_vol = area * dz
+
+# create a new reactor
+r2 = ct.IdealGasReactor(gas2)
+r2.volume = r_vol
+
+# create a reservoir to represent the reactor immediately upstream. Note
+# that the gas object is set already to the state of the upstream reactor
+upstream = ct.Reservoir(gas2, name='upstream')
+
+# create a reservoir for the reactor to exhaust into. The composition of
+# this reservoir is irrelevant.
+downstream = ct.Reservoir(gas2, name='downstream')
+
+# The mass flow rate into the reactor will be fixed by using a
+# MassFlowController object.
+m = ct.MassFlowController(upstream, r2, mdot=mass_flow_rate2)
+
+# We need an outlet to the downstream reservoir. This will determine the
+# pressure in the reactor. The value of K will only affect the transient
+# pressure difference.
+v = ct.PressureController(r2, downstream, master=m, K=1e-5)
+
+sim2 = ct.ReactorNet([r2])
+
+# define time, space, and other information vectors
+z2 = (np.arange(n_steps) + 1) * dz
+t_r2 = np.zeros_like(z2)  # residence time in each reactor
+t2 = np.zeros_like(z2)
+u2 = np.zeros_like(z2)
+T2 = np.zeros_like(z2)
+X_H2_2 = np.zeros_like(z2)
+# iterate through the PFR cells
+for n in range(n_steps):
+    # Set the state of the reservoir to match that of the previous reactor
+    gas2.TDY = r2.thermo.TDY
+    upstream.syncState()
+    # integrate the reactor forward in time until steady state is reached
+    sim2.set_initial_time(0)  # forces reinitialization
+    time = 0
+    all_done = False
+    # determine steady state from H2 mole fraction
+    X_H2_previous = r2.thermo['H2'].X
+    while not all_done:
+        time += dt
+        sim2.advance(time)
+        if np.abs(r2.thermo['H2'].X - X_H2_previous) < 1.e-10:
+            # check whether surface coverages are in steady state.
+            all_done = True
+        else:
+            X_H2_previous = r2.thermo['H2'].X
+    # compute velocity and transform into time
+    u2[n] = mass_flow_rate2 / area / r2.thermo.density
+    t_r2[n] = r2.mass / mass_flow_rate2  # residence time in this reactor
+    t2[n] = np.sum(t_r2)
+    # write output data
+    T2[n] = r2.T
+    X_H2_2[n] = r2.thermo['H2'].X
+
+#####################################################################
+
+
+#####################################################################
+# Compare Results in matplotlib
+#####################################################################
+
+import matplotlib.pyplot as plt
+
+plt.figure()
+plt.plot(z1, T1, label='Lagrangian Particle')
+plt.plot(z2, T2, label='Reactor Chain')
+plt.xlabel('$z$ [m]')
+plt.ylabel('$T$ [K]')
+plt.legend(loc=0)
+plt.show()
+plt.savefig('pfr_T_z.png')
+
+plt.figure()
+plt.plot(t1, X_H2_1, label='Lagrangian Particle')
+plt.plot(t2, X_H2_2, label='Reactor Chain')
+plt.xlabel('$t$ [s]')
+plt.ylabel('$X_{H_2}$ [-]')
+plt.legend(loc=0)
+plt.show()
+plt.savefig('pfr_XH2_t.png')