diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C
index c5dec29661..e5d3026828 100644
--- a/applications/test/ODE/Test-ODE.C
+++ b/applications/test/ODE/Test-ODE.C
@@ -27,7 +27,7 @@ Description
 
 #include "argList.H"
 #include "IOmanip.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "ODESolver.H"
 #include "RK.H"
 
@@ -37,7 +37,7 @@ using namespace Foam;
 
 class testODE
 :
-    public ODE
+    public ODESystem
 {
 
 public:
@@ -107,9 +107,13 @@ int main(int argc, char *argv[])
     argList::validArgs.append("ODESolver");
     argList args(argc, argv);
 
+    // Create the ODE system
     testODE ode;
+
+    // Create the selected ODE system solver
     autoPtr<ODESolver> odeSolver = ODESolver::New(args[1], ode);
 
+    // Initialise the ODE system fields
     scalar xStart = 1.0;
     scalarField yStart(ode.nEqns());
     yStart[0] = ::Foam::j0(xStart);
@@ -117,6 +121,7 @@ int main(int argc, char *argv[])
     yStart[2] = ::Foam::jn(2, xStart);
     yStart[3] = ::Foam::jn(3, xStart);
 
+    // Print the evolution of the solution and the time-step
     scalarField dyStart(ode.nEqns());
     ode.derivatives(xStart, yStart, dyStart);
 
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C
index 9246605b51..e9218641a7 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.C
+++ b/src/ODE/ODESolvers/KRR4/KRR4.C
@@ -55,7 +55,7 @@ const scalar
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::KRR4::KRR4(const ODE& ode)
+Foam::KRR4::KRR4(const ODESystem& ode)
 :
     ODESolver(ode),
     yTemp_(n_, 0.0),
@@ -76,7 +76,7 @@ Foam::KRR4::KRR4(const ODE& ode)
 
 void Foam::KRR4::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
@@ -168,7 +168,7 @@ void Foam::KRR4::solve
             (
                 "void Foam::KRR4::solve"
                 "("
-                    "const ODE&, "
+                    "const ODESystem&, "
                     "scalar&, "
                     "scalarField&, "
                     "scalarField&, "
@@ -206,7 +206,7 @@ void Foam::KRR4::solve
     (
         "void Foam::KRR4::solve"
         "("
-            "const ODE&, "
+            "const ODESystem&, "
             "scalar&, "
             "scalarField&, "
             "scalarField&, "
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.H b/src/ODE/ODESolvers/KRR4/KRR4.H
index 17e9a8caa0..1eb434a0ba 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.H
+++ b/src/ODE/ODESolvers/KRR4/KRR4.H
@@ -88,14 +88,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        KRR4(const ODE& ode);
+        KRR4(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C
index 7a981dfd16..5f7188a079 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.C
+++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C
@@ -36,7 +36,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::ODESolver::ODESolver(const ODE& ode)
+Foam::ODESolver::ODESolver(const ODESystem& ode)
 :
     n_(ode.nEqns()),
     yScale_(n_),
@@ -48,7 +48,7 @@ Foam::ODESolver::ODESolver(const ODE& ode)
 
 void Foam::ODESolver::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
@@ -102,7 +102,7 @@ void Foam::ODESolver::solve
     FatalErrorIn
     (
         "ODESolver::solve"
-        "(const ODE& ode, const scalar xStart, const scalar xEnd,"
+        "(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
         "scalarField& yStart, const scalar eps, scalar& hEst) const"
     )   << "Too many integration steps"
         << exit(FatalError);
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H
index 79c3fc5d1b..04315ea7dc 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.H
+++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::ODESolver
 
 Description
-    Selection for ODE solver
+    Abstract base-class for ODE system solvers
 
 SourceFiles
     ODESolver.C
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef ODESolver_H
 #define ODESolver_H
 
-#include "ODE.H"
+#include "ODESystem.H"
 #include "typeInfo.H"
 #include "autoPtr.H"
 
@@ -82,7 +82,7 @@ public:
             autoPtr,
             ODESolver,
             ODE,
-            (const ODE& ode),
+            (const ODESystem& ode),
             (ode)
         );
 
@@ -90,7 +90,7 @@ public:
     // Constructors
 
         //- Construct for given ODE
-        ODESolver(const ODE& ode);
+        ODESolver(const ODESystem& ode);
 
 
     // Selectors
@@ -99,7 +99,7 @@ public:
         static autoPtr<ODESolver> New
         (
             const word& ODESolverTypeName,
-            const ODE& ode
+            const ODESystem& ode
         );
 
 
@@ -112,7 +112,7 @@ public:
 
         virtual void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
@@ -126,7 +126,7 @@ public:
 
         virtual void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar xStart,
             const scalar xEnd,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C b/src/ODE/ODESolvers/ODESolver/ODESolverNew.C
index 186e557726..ab4c3a0cb7 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C
+++ b/src/ODE/ODESolvers/ODESolver/ODESolverNew.C
@@ -30,7 +30,7 @@ License
 Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
 (
     const Foam::word& ODESolverTypeName,
-    const Foam::ODE& ode
+    const Foam::ODESystem& ode
 )
 {
     Info<< "Selecting ODE solver " << ODESolverTypeName << endl;
@@ -42,7 +42,7 @@ Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
     {
         FatalErrorIn
         (
-            "ODESolver::New(const word& ODESolverTypeName, const ODE& ode)"
+            "ODESolver::New(const word& ODESolverTypeName, const ODESystem& ode)"
         )   << "Unknown ODESolver type "
             << ODESolverTypeName << nl << nl
             << "Valid ODESolvers are : " << endl
diff --git a/src/ODE/ODESolvers/RK/RK.C b/src/ODE/ODESolvers/RK/RK.C
index 6ae5036078..06c5e6ae81 100644
--- a/src/ODE/ODESolvers/RK/RK.C
+++ b/src/ODE/ODESolvers/RK/RK.C
@@ -54,7 +54,7 @@ const scalar
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::RK::RK(const ODE& ode)
+Foam::RK::RK(const ODESystem& ode)
 :
     ODESolver(ode),
     yTemp_(n_, 0.0),
@@ -72,7 +72,7 @@ Foam::RK::RK(const ODE& ode)
 
 void Foam::RK::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar x,
     const scalarField& y,
     const scalarField& dydx,
@@ -142,7 +142,7 @@ void Foam::RK::solve
 
 void Foam::RK::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
diff --git a/src/ODE/ODESolvers/RK/RK.H b/src/ODE/ODESolvers/RK/RK.H
index 5ab15fa0d1..20fa3aa789 100644
--- a/src/ODE/ODESolvers/RK/RK.H
+++ b/src/ODE/ODESolvers/RK/RK.H
@@ -81,14 +81,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        RK(const ODE& ode);
+        RK(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar x,
             const scalarField& y,
             const scalarField& dydx,
@@ -100,7 +100,7 @@ public:
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C
index d86f60d5c9..95ab44776c 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.C
+++ b/src/ODE/ODESolvers/SIBS/SIBS.C
@@ -47,7 +47,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::SIBS::SIBS(const ODE& ode)
+Foam::SIBS::SIBS(const ODESystem& ode)
 :
     ODESolver(ode),
     a_(iMaxX_, 0.0),
@@ -70,7 +70,7 @@ Foam::SIBS::SIBS(const ODE& ode)
 
 void Foam::SIBS::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H
index 7c1d9c6fbe..a031b1c83b 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.H
+++ b/src/ODE/ODESolvers/SIBS/SIBS.H
@@ -78,7 +78,7 @@ class SIBS
 
         void SIMPR
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar xStart,
             const scalarField& y,
             const scalarField& dydx,
@@ -110,14 +110,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        SIBS(const ODE& ode);
+        SIBS(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C
index 9c67ce1158..36402a202d 100644
--- a/src/ODE/ODESolvers/SIBS/SIMPR.C
+++ b/src/ODE/ODESolvers/SIBS/SIMPR.C
@@ -29,7 +29,7 @@ License
 
 void Foam::SIBS::SIMPR
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar xStart,
     const scalarField& y,
     const scalarField& dydx,
diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODESystem/ODESystem.H
similarity index 80%
rename from src/ODE/ODE/ODE.H
rename to src/ODE/ODESystem/ODESystem.H
index 4c98d8dc4c..90da37c427 100644
--- a/src/ODE/ODE/ODE.H
+++ b/src/ODE/ODESystem/ODESystem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,15 +22,15 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::ODE
+    Foam::ODESystem
 
 Description
-    Abstract base class for the ODE solvers.
+    Abstract base class for the systems of ordinary differential equations.
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ODE_H
-#define ODE_H
+#ifndef ODESystem_H
+#define ODESystem_H
 
 #include "scalarField.H"
 #include "scalarMatrices.H"
@@ -41,30 +41,32 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class ODE Declaration
+                       Class ODESystem Declaration
 \*---------------------------------------------------------------------------*/
 
-class ODE
+class ODESystem
 {
 
 public:
 
     // Constructors
 
-        //- Construct null
-        ODE()
+        //- Construct null8
+        ODESystem()
         {}
 
 
     //- Destructor
-    virtual ~ODE()
+    virtual ~ODESystem()
     {}
 
 
     // Member Functions
 
+        //- Return the number of equations in the system
         virtual label nEqns() const = 0;
 
+        //- Calculate the derivatives in dydx
         virtual void derivatives
         (
             const scalar x,
@@ -72,7 +74,8 @@ public:
             scalarField& dydx
         ) const = 0;
 
-
+        //- Calculate the Jacobian of the system
+        //  Need by the stiff-system solvers
         virtual void jacobian
         (
             const scalar x,
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index 6d5ad3d562..8357df6dff 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -36,7 +36,7 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
 )
 :
     CompType(mesh),
-    ODE(),
+    ODESystem(),
     Y_(this->thermo().composition().Y()),
     reactions_
     (
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
index 6da538d3d9..87bba8c2bb 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@@ -39,7 +39,7 @@ SourceFiles
 #define chemistryModel_H
 
 #include "Reaction.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "volFieldsFwd.H"
 #include "simpleMatrix.H"
 #include "DimensionedField.H"
@@ -60,7 +60,7 @@ template<class CompType, class ThermoType>
 class chemistryModel
 :
     public CompType,
-    public ODE
+    public ODESystem
 {
     // Private Member Functions
 
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
index b60b5c7eb5..802a8ab87d 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
@@ -37,7 +37,7 @@ solidChemistryModel
 )
 :
     CompType(mesh),
-    ODE(),
+    ODESystem(),
     Ys_(this->solidThermo().composition().Y()),
     reactions_
     (
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
index 09b04d06a9..23cbbde6b7 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
@@ -40,7 +40,7 @@ SourceFiles
 #define solidChemistryModel_H
 
 #include "Reaction.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "volFieldsFwd.H"
 #include "DimensionedField.H"
 #include "simpleMatrix.H"
@@ -61,7 +61,7 @@ template<class CompType, class SolidThermo>
 class solidChemistryModel
 :
     public CompType,
-    public ODE
+    public ODESystem
 {
     // Private Member Functions