diff --git a/NumLib/ODESolver/TimeDiscretization.cpp b/NumLib/ODESolver/TimeDiscretization.cpp
index 820d9204bf0bd47630b1a36ca3af7ef2e56689cc..99511caa4d368aacbafc4ad3c351cb78fe80d9c3 100644
--- a/NumLib/ODESolver/TimeDiscretization.cpp
+++ b/NumLib/ODESolver/TimeDiscretization.cpp
@@ -47,4 +47,28 @@ double computeRelativeChangeFromPreviousTimestep(GlobalVector const& x,
     // Only norm_x is close to zero
     return norm_dx / std::numeric_limits<double>::epsilon();
 }
+
+void TimeDiscretization::getXdot(GlobalVector const& x_at_new_timestep,
+                                 GlobalVector const& x_old,
+                                 GlobalVector& xdot) const
+{
+    namespace LinAlg = MathLib::LinAlg;
+
+    double const dt = getCurrentTimeIncrement();
+
+    // xdot = 1/dt * x_at_new_timestep - x_old
+    getWeightedOldX(xdot, x_old);
+    LinAlg::axpby(xdot, 1. / dt, -1.0, x_at_new_timestep);
+}
+
+void BackwardEuler::getWeightedOldX(GlobalVector& y,
+                                    GlobalVector const& x_old) const
+{
+    namespace LinAlg = MathLib::LinAlg;
+
+    // y = x_old / delta_t
+    LinAlg::copy(x_old, y);
+    LinAlg::scale(y, 1.0 / _delta_t);
+}
+
 }  // end of namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index 2ae33459c1e16af192235cac5ef438887e76aa71..6ad98e921ce518263cf83174710fabab6f96795b 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -132,16 +132,7 @@ public:
     //! \f$.
     void getXdot(GlobalVector const& x_at_new_timestep,
                  GlobalVector const& x_old,
-                 GlobalVector& xdot) const
-    {
-        namespace LinAlg = MathLib::LinAlg;
-
-        double const dt = getCurrentTimeIncrement();
-
-        // xdot = 1/dt * x_at_new_timestep - x_old
-        getWeightedOldX(xdot, x_old);
-        LinAlg::axpby(xdot, 1. / dt, -1.0, x_at_new_timestep);
-    }
+                 GlobalVector& xdot) const;
 
     //! Returns \f$ x_O \f$.
     virtual void getWeightedOldX(
@@ -155,7 +146,6 @@ class BackwardEuler final : public TimeDiscretization
 {
 public:
     void setInitialState(const double t0) override { _t = t0; }
-
     void nextTimestep(const double t, const double delta_t) override
     {
         _t = t;
@@ -165,14 +155,7 @@ public:
     double getCurrentTime() const override { return _t; }
     double getCurrentTimeIncrement() const override { return _delta_t; }
     void getWeightedOldX(GlobalVector& y,
-                         GlobalVector const& x_old) const override
-    {
-        namespace LinAlg = MathLib::LinAlg;
-
-        // y = x_old / delta_t
-        LinAlg::copy(x_old, y);
-        LinAlg::scale(y, 1.0 / _delta_t);
-    }
+                         GlobalVector const& x_old) const override;
 
 private:
     double _t = std::numeric_limits<double>::quiet_NaN();  //!< \f$ t_C \f$