diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C
index 4222fa4cc4..ae20374259 100644
--- a/applications/test/ODE/Test-ODE.C
+++ b/applications/test/ODE/Test-ODE.C
@@ -142,7 +142,7 @@ int main(int argc, char *argv[])
         scalar dxNext = dxEst;
 
         odeSolver->relTol() = relTol;
-        odeSolver->solve(ode, x, y, dxNext);
+        odeSolver->solve(x, y, dxNext);
 
         Info<< scientific << setw(13) << relTol;
         Info<< fixed << setw(11) << dxEst;
@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
     scalar dxEst = 0.5;
 
     odeSolver->relTol() = 1e-4;
-    odeSolver->solve(ode, x, xEnd, y, dxEst);
+    odeSolver->solve(x, xEnd, y, dxEst);
 
     Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
     Info      << "Numerical:  y(2.0) = " << y << ", dxEst = " << dxEst << endl;
diff --git a/src/ODE/ODESolvers/Euler/Euler.C b/src/ODE/ODESolvers/Euler/Euler.C
index 1a6edf576c..1496099ae8 100644
--- a/src/ODE/ODESolvers/Euler/Euler.C
+++ b/src/ODE/ODESolvers/Euler/Euler.C
@@ -49,7 +49,6 @@ Foam::Euler::Euler(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Euler::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -75,13 +74,12 @@ Foam::scalar Foam::Euler::solve
 
 void Foam::Euler::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/Euler/Euler.H b/src/ODE/ODESolvers/Euler/Euler.H
index 90f67f4b4c..3cefc425be 100644
--- a/src/ODE/ODESolvers/Euler/Euler.H
+++ b/src/ODE/ODESolvers/Euler/Euler.H
@@ -84,7 +84,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -95,7 +94,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.C b/src/ODE/ODESolvers/EulerSI/EulerSI.C
index 1785101758..9b28f733ea 100644
--- a/src/ODE/ODESolvers/EulerSI/EulerSI.C
+++ b/src/ODE/ODESolvers/EulerSI/EulerSI.C
@@ -54,7 +54,6 @@ Foam::EulerSI::EulerSI(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::EulerSI::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -62,7 +61,7 @@ Foam::scalar Foam::EulerSI::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -95,13 +94,12 @@ Foam::scalar Foam::EulerSI::solve
 
 void Foam::EulerSI::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.H b/src/ODE/ODESolvers/EulerSI/EulerSI.H
index 392e818e87..656e3df0a4 100644
--- a/src/ODE/ODESolvers/EulerSI/EulerSI.H
+++ b/src/ODE/ODESolvers/EulerSI/EulerSI.H
@@ -91,7 +91,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -102,7 +101,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C
index f4dd3588eb..3f7f0e0996 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.C
+++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C
@@ -38,6 +38,7 @@ namespace Foam
 
 Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict)
 :
+    odes_(ode),
     n_(ode.nEqns()),
     absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
     relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-4)),
@@ -52,6 +53,7 @@ Foam::ODESolver::ODESolver
     const scalarField& relTol
 )
 :
+    odes_(ode),
     n_(ode.nEqns()),
     absTol_(absTol),
     relTol_(relTol),
@@ -82,7 +84,6 @@ Foam::scalar Foam::ODESolver::normalizeError
 
 void Foam::ODESolver::solve
 (
-    const ODESystem& ode,
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
@@ -105,7 +106,7 @@ void Foam::ODESolver::solve
         }
 
         // Integrate as far as possible up to dxEst
-        solve(ode, x, y, dxEst);
+        solve(x, y, dxEst);
 
         // Check if reached xEnd
         if ((x - xEnd)*(xEnd - xStart) >= 0)
@@ -122,7 +123,7 @@ void Foam::ODESolver::solve
     FatalErrorIn
     (
         "ODESolver::solve"
-        "(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
+        "(const scalar xStart, const scalar xEnd,"
         "scalarField& y, scalar& dxEst) const"
     )   << "Integration steps greater than maximum " << maxSteps_
         << exit(FatalError);
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H
index bd767a9063..4c3becd507 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.H
+++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H
@@ -55,6 +55,9 @@ protected:
 
     // Protected data
 
+        //- Reference to ODESystem
+        const ODESystem& odes_;
+
         //- Size of the ODESystem
         label n_;
 
@@ -150,7 +153,6 @@ public:
         //  Update the state and return an estimate for the next step in dxTry
         virtual void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
@@ -160,7 +162,6 @@ public:
         //  and return an estimate for the next step in dxTry
         virtual void solve
         (
-            const ODESystem& ode,
             const scalar xStart,
             const scalar xEnd,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.C b/src/ODE/ODESolvers/RKCK45/RKCK45.C
index 2de3cdf99e..30ae91de1a 100644
--- a/src/ODE/ODESolvers/RKCK45/RKCK45.C
+++ b/src/ODE/ODESolvers/RKCK45/RKCK45.C
@@ -89,7 +89,6 @@ Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKCK45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -102,21 +101,21 @@ Foam::scalar Foam::RKCK45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -124,7 +123,7 @@ Foam::scalar Foam::RKCK45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -133,7 +132,7 @@ Foam::scalar Foam::RKCK45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + c6*dx, yTemp_, k6_);
+    odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
 
     forAll(y, i)
     {
@@ -154,13 +153,12 @@ Foam::scalar Foam::RKCK45::solve
 
 void Foam::RKCK45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.H b/src/ODE/ODESolvers/RKCK45/RKCK45.H
index a6b551a475..12edcf982a 100644
--- a/src/ODE/ODESolvers/RKCK45/RKCK45.H
+++ b/src/ODE/ODESolvers/RKCK45/RKCK45.H
@@ -110,7 +110,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -121,7 +120,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.C b/src/ODE/ODESolvers/RKDP45/RKDP45.C
index cf319dca8b..f9a3a78836 100644
--- a/src/ODE/ODESolvers/RKDP45/RKDP45.C
+++ b/src/ODE/ODESolvers/RKDP45/RKDP45.C
@@ -94,7 +94,6 @@ Foam::RKDP45::RKDP45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKDP45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -107,21 +106,21 @@ Foam::scalar Foam::RKDP45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -129,7 +128,7 @@ Foam::scalar Foam::RKDP45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -138,7 +137,7 @@ Foam::scalar Foam::RKDP45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + dx, yTemp_, k6_);
+    odes_.derivatives(x0 + dx, yTemp_, k6_);
 
     forAll(y, i)
     {
@@ -147,7 +146,7 @@ Foam::scalar Foam::RKDP45::solve
     }
 
     // Reuse k2_ for the derivative of the new state
-    odes.derivatives(x0 + dx, y, k2_);
+    odes_.derivatives(x0 + dx, y, k2_);
 
     forAll(err_, i)
     {
@@ -163,13 +162,12 @@ Foam::scalar Foam::RKDP45::solve
 
 void Foam::RKDP45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.H b/src/ODE/ODESolvers/RKDP45/RKDP45.H
index 95ae212e66..40e23fafe8 100644
--- a/src/ODE/ODESolvers/RKDP45/RKDP45.H
+++ b/src/ODE/ODESolvers/RKDP45/RKDP45.H
@@ -114,7 +114,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -125,7 +124,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/RKF45/RKF45.C b/src/ODE/ODESolvers/RKF45/RKF45.C
index 661c618255..92804e0634 100644
--- a/src/ODE/ODESolvers/RKF45/RKF45.C
+++ b/src/ODE/ODESolvers/RKF45/RKF45.C
@@ -90,7 +90,6 @@ Foam::RKF45::RKF45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKF45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -103,21 +102,21 @@ Foam::scalar Foam::RKF45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -125,7 +124,7 @@ Foam::scalar Foam::RKF45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -134,7 +133,7 @@ Foam::scalar Foam::RKF45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + c6*dx, yTemp_, k6_);
+    odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
 
     // Calculate the 5th-order solution
     forAll(y, i)
@@ -159,13 +158,12 @@ Foam::scalar Foam::RKF45::solve
 
 void Foam::RKF45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/RKF45/RKF45.H b/src/ODE/ODESolvers/RKF45/RKF45.H
index df91ed7808..1f6cfd1fe9 100644
--- a/src/ODE/ODESolvers/RKF45/RKF45.H
+++ b/src/ODE/ODESolvers/RKF45/RKF45.H
@@ -138,7 +138,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -149,7 +148,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C
index 13371a8e90..15e81ed3ea 100644
--- a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C
+++ b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C
@@ -68,7 +68,6 @@ Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Rosenbrock21::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -76,7 +75,7 @@ Foam::scalar Foam::Rosenbrock21::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -104,7 +103,7 @@ Foam::scalar Foam::Rosenbrock21::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -126,13 +125,12 @@ Foam::scalar Foam::Rosenbrock21::solve
 
 void Foam::Rosenbrock21::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H
index a29207e129..3c3c40825d 100644
--- a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H
+++ b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H
@@ -101,7 +101,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -112,7 +111,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C
index e471e3ffe1..0eba13b97c 100644
--- a/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C
+++ b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C
@@ -81,7 +81,6 @@ Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Rosenbrock32::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -89,7 +88,7 @@ Foam::scalar Foam::Rosenbrock32::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -117,7 +116,7 @@ Foam::scalar Foam::Rosenbrock32::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -148,13 +147,12 @@ Foam::scalar Foam::Rosenbrock32::solve
 
 void Foam::Rosenbrock32::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H
index 8cd2e50dbd..596ab2ad3e 100644
--- a/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H
+++ b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H
@@ -105,7 +105,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -116,7 +115,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C
index 061390fdbc..511261792a 100644
--- a/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C
+++ b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C
@@ -126,7 +126,6 @@ Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Rosenbrock43::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -134,7 +133,7 @@ Foam::scalar Foam::Rosenbrock43::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -162,7 +161,7 @@ Foam::scalar Foam::Rosenbrock43::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -177,7 +176,7 @@ Foam::scalar Foam::Rosenbrock43::solve
         y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
     }
 
-    ode.derivatives(x0 + c3*dx, y, dydx_);
+    odes_.derivatives(x0 + c3*dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -208,13 +207,12 @@ Foam::scalar Foam::Rosenbrock43::solve
 
 void Foam::Rosenbrock43::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.H b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.H
index f519f84b35..ef84759c37 100644
--- a/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.H
+++ b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.H
@@ -112,7 +112,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -123,7 +122,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C
index eae8552a47..b6b6cf2c1d 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.C
+++ b/src/ODE/ODESolvers/SIBS/SIBS.C
@@ -70,13 +70,12 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
 
 void Foam::SIBS::solve
 (
-    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    ode.derivatives(x, y, dydx0_);
+    odes_.derivatives(x, y, dydx0_);
 
     scalar h = dxTry;
     bool exitflag = false;
@@ -124,7 +123,7 @@ void Foam::SIBS::solve
     label k = 0;
     yTemp_ = y;
 
-    ode.jacobian(x, y, dfdx_, dfdy_);
+    odes_.jacobian(x, y, dfdx_, dfdy_);
 
     if (x != xNew_ || h != dxTry)
     {
@@ -151,7 +150,7 @@ void Foam::SIBS::solve
                     << exit(FatalError);
             }
 
-            SIMPR(ode, x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
+            SIMPR(x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
             scalar xest = sqr(h/nSeq_[k]);
 
             polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H
index 03664fd53a..0ef0a1d40c 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.H
+++ b/src/ODE/ODESolvers/SIBS/SIBS.H
@@ -84,7 +84,6 @@ class SIBS
 
         void SIMPR
         (
-            const ODESystem& ode,
             const scalar xStart,
             const scalarField& y,
             const scalarField& dydx,
@@ -123,7 +122,6 @@ public:
 
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C
index 36402a202d..d21a931e1d 100644
--- a/src/ODE/ODESolvers/SIBS/SIMPR.C
+++ b/src/ODE/ODESolvers/SIBS/SIMPR.C
@@ -29,7 +29,6 @@ License
 
 void Foam::SIBS::SIMPR
 (
-    const ODESystem& ode,
     const scalar xStart,
     const scalarField& y,
     const scalarField& dydx,
@@ -72,7 +71,7 @@ void Foam::SIBS::SIMPR
 
     scalar x = xStart + h;
 
-    ode.derivatives(x, ytemp, yEnd);
+    odes_.derivatives(x, ytemp, yEnd);
 
     for (register label nn=2; nn<=nSteps; nn++)
     {
@@ -90,7 +89,7 @@ void Foam::SIBS::SIMPR
 
         x += h;
 
-        ode.derivatives(x, ytemp, yEnd);
+        odes_.derivatives(x, ytemp, yEnd);
     }
     for (register label i=0; i<n_; i++)
     {
diff --git a/src/ODE/ODESolvers/Trapezoid/Trapezoid.C b/src/ODE/ODESolvers/Trapezoid/Trapezoid.C
index 8a2c281433..027308434d 100644
--- a/src/ODE/ODESolvers/Trapezoid/Trapezoid.C
+++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.C
@@ -49,7 +49,6 @@ Foam::Trapezoid::Trapezoid(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Trapezoid::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -64,7 +63,7 @@ Foam::scalar Foam::Trapezoid::solve
     }
 
     // Evaluate the system for the predicted state
-    ode.derivatives(x0 + dx, y, err_);
+    odes_.derivatives(x0 + dx, y, err_);
 
     // Update the state as the average between the prediction and the correction
     // and estimate the error from the difference
@@ -80,13 +79,12 @@ Foam::scalar Foam::Trapezoid::solve
 
 void Foam::Trapezoid::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/Trapezoid/Trapezoid.H b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H
index 327dbe44bc..c52d80f51d 100644
--- a/src/ODE/ODESolvers/Trapezoid/Trapezoid.H
+++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H
@@ -74,7 +74,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -85,7 +84,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C
index 13951ae7b3..ee49a70d88 100644
--- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C
+++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C
@@ -64,7 +64,7 @@ void Foam::adaptiveSolver::solve
     do
     {
         // Solve step and provide error estimate
-        err = solve(odes, x, y, dydx0_, dx, yTemp_);
+        err = solve(x, y, dydx0_, dx, yTemp_);
 
         // If error is large reduce dx
         if (err > 1)
diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H
index d86bb88ce1..0d7dcb035c 100644
--- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H
+++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H
@@ -77,7 +77,6 @@ public:
         //- Solve a single step dx and return the error
         virtual scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
diff --git a/src/ODE/ODESolvers/rodas32/rodas32.C b/src/ODE/ODESolvers/rodas32/rodas32.C
index 6221788ea6..fb4479b0b0 100644
--- a/src/ODE/ODESolvers/rodas32/rodas32.C
+++ b/src/ODE/ODESolvers/rodas32/rodas32.C
@@ -72,7 +72,6 @@ Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::rodas32::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -80,7 +79,7 @@ Foam::scalar Foam::rodas32::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -117,7 +116,7 @@ Foam::scalar Foam::rodas32::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -133,7 +132,7 @@ Foam::scalar Foam::rodas32::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(err_, i)
     {
@@ -153,13 +152,12 @@ Foam::scalar Foam::rodas32::solve
 
 void Foam::rodas32::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/rodas32/rodas32.H b/src/ODE/ODESolvers/rodas32/rodas32.H
index 5cd05e9a65..f3628be612 100644
--- a/src/ODE/ODESolvers/rodas32/rodas32.H
+++ b/src/ODE/ODESolvers/rodas32/rodas32.H
@@ -105,7 +105,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -116,7 +115,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
diff --git a/src/ODE/ODESolvers/rodas43/rodas43.C b/src/ODE/ODESolvers/rodas43/rodas43.C
index 92b78a7ab4..3caf0854ce 100644
--- a/src/ODE/ODESolvers/rodas43/rodas43.C
+++ b/src/ODE/ODESolvers/rodas43/rodas43.C
@@ -95,7 +95,6 @@ Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::rodas43::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -103,7 +102,7 @@ Foam::scalar Foam::rodas43::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -131,7 +130,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -146,7 +145,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
     }
 
-    ode.derivatives(x0 + c3*dx, y, dydx_);
+    odes_.derivatives(x0 + c3*dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -161,7 +160,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
     }
 
-    ode.derivatives(x0 + c4*dx, y, dydx_);
+    odes_.derivatives(x0 + c4*dx, y, dydx_);
 
     forAll(k4_, i)
     {
@@ -178,7 +177,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(k5_, i)
     {
@@ -195,7 +194,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(err_, i)
     {
@@ -216,13 +215,12 @@ Foam::scalar Foam::rodas43::solve
 
 void Foam::rodas43::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 
diff --git a/src/ODE/ODESolvers/rodas43/rodas43.H b/src/ODE/ODESolvers/rodas43/rodas43.H
index e77c0b423b..21bb3d5b66 100644
--- a/src/ODE/ODESolvers/rodas43/rodas43.H
+++ b/src/ODE/ODESolvers/rodas43/rodas43.H
@@ -106,7 +106,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -117,7 +116,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry