diff --git a/src/cellmodels/MBeelerReuter.cpp b/src/cellmodels/MBeelerReuter.cpp
index 3d243fc095043fa1ea5a9083b3d7c0b8bcf98796..0323c9033abb2ae2dee8cfc1d962faf286985a80 100644
--- a/src/cellmodels/MBeelerReuter.cpp
+++ b/src/cellmodels/MBeelerReuter.cpp
@@ -15,38 +15,6 @@
 #define BR_yinfty_m(V)   (Alpha(4,V) * BR_tau_m(V)  )
 #define BR_yinfty_x_1(V) (Alpha(5,V) * BR_tau_x_1(V))
 
-//w=(d,f,h,j,m,x_1)
-//c=(Ca)
-void
-MBeelerReuter::EvaluateRHS(double t, const std::vector<double> &var, std::vector<double> &g) const {
-
-  double V = BR_Pot(var);
-  double I_ext = IExt(t);
-
-
-  double I_Na =
-      (g_Na * BR_m(var) * BR_m(var) * BR_m(var) * BR_h(var) * BR_j(var) + g_NaC) * (V - E_Na);
-  double I_s = g_s * BR_d(var) * BR_f(var) *
-               (V - (-82.3 - 13.0287 * log(BR_Ca(var))));//Achtung, ist log=ln????
-  double I_x1 = BR_x_1(var) * ((0.8 * (exp(0.04 * (V + 77)) - 1)) / (exp(0.04 * (V + 35))));
-  double I_K = 0.35 * (((4 * (exp(0.04 * (V + 85)) - 1)) /
-                        (exp(0.08 * (V + 53)) + exp(0.04 * (V + 53)))) +
-                       ((0.2 * (V + 23)) / (1 - exp(-0.04 * (V + 23)))));
-
-
-  BR_Pot(g) = -1.0 * (I_Na + I_s + I_x1 + I_K - I_ext);//*0.001
-
-  BR_Ca(g) = -0.0000001 * I_s + 0.07 * (0.0000001 - BR_Ca(var));
-
-  BR_d(g) = (BR_yinfty_d(V) - BR_d(var)) / BR_tau_d(V);
-  BR_f(g) = (BR_yinfty_f(V) - BR_f(var)) / BR_tau_f(V);
-  BR_h(g) = (BR_yinfty_h(V) - BR_h(var)) / BR_tau_h(V);
-  BR_j(g) = (BR_yinfty_j(V) - BR_j(var)) / BR_tau_j(V);
-  BR_m(g) = (BR_yinfty_m(V) - BR_m(var)) / BR_tau_m(V);
-  BR_x_1(g) = (BR_yinfty_x_1(V) - BR_x_1(var)) / BR_tau_x_1(V);
-
-}
-
 double MBeelerReuter::GetIion(const std::vector<double> &var) const {
   double mVolt = var[0];
   double Ca = var[1];
diff --git a/src/cellmodels/MBeelerReuter.hpp b/src/cellmodels/MBeelerReuter.hpp
index 90ea63c3755a0ba370fcc0b673df5fb28bd7712a..7b540d4b8c80643c481f7bc4ffac72a8abd9163a 100644
--- a/src/cellmodels/MBeelerReuter.hpp
+++ b/src/cellmodels/MBeelerReuter.hpp
@@ -97,9 +97,6 @@ public:
   double GetIion(const std::vector<double> &var) const override;
 
   double JacobiEntry(int i, int j, const std::vector<double> &VW) const override;
-
-  void
-  EvaluateRHS(double t, const std::vector<double> &var, std::vector<double> &g) const override;
 };
 
 #endif //MBEELERREUTER_HPP
diff --git a/src/cellmodels/MCellModel.cpp b/src/cellmodels/MCellModel.cpp
index b4af2b4c14be8c8e516e1d836ba04d8cadab97aa..76050a27ac413e2263379cec9f8397e56a12e061 100644
--- a/src/cellmodels/MCellModel.cpp
+++ b/src/cellmodels/MCellModel.cpp
@@ -76,10 +76,15 @@ void MTensionModel::EvaluateRHS(double t, const vector<double> &var, vector<doub
 }*/
 
 
-void MCellModel::RHSUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
+void MElphyModel::RHSUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
   IIonUpdate(var, g, iExt);
   ConcentrationUpdate(var, g, iExt);
   GatingUpdate(var, g);
 }
 
+void MTensionModel::RHSUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
+  elphy->RHSUpdate(var, g, iExt);
+  TensionUpdate(var, g, iExt);
+}
+
 
diff --git a/src/cellmodels/MCellModel.hpp b/src/cellmodels/MCellModel.hpp
index 06efc46f9102c345f232d3bae36045c6caff17f1..b1354ae5391e43bb702d5b72b5962c6f678f31d7 100644
--- a/src/cellmodels/MCellModel.hpp
+++ b/src/cellmodels/MCellModel.hpp
@@ -37,7 +37,7 @@ public:
 
   virtual int CalciumIndex() const { return 0; }
 
-  double IExt(double t) const {
+  virtual double IExt(double t) const {
     return excitationFunc(t, amplitude, excitation, duration);
   }
 
@@ -61,7 +61,7 @@ public:
   ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const {};
 
   virtual void
-  RHSUpdate(const std::vector<double> &var, vector<double> &g, double iExt) const;
+  RHSUpdate(const std::vector<double> &var, vector<double> &g, double iExt) const = 0;
 
   virtual void
   IIonUpdate(const std::vector<double> &var, vector<double> &g, double iExt) const {};
@@ -71,10 +71,7 @@ public:
 
   virtual void GatingUpdate(const std::vector<double> &var, vector<double> &g) const {};
 
-  virtual void UpdateValues(double dt, std::vector<double> &var, const vector<double> &g) const {};
-
-  virtual void
-  EvaluateRHS(double t, const std::vector<double> &var, std::vector<double> &g) const = 0;
+  virtual void UpdateValues(double dt, std::vector<double> &var, const vector<double> &g) const = 0;
 
   virtual double InitialValue() const {Exit("Not implemented in the basis class")};
 
@@ -137,6 +134,7 @@ public:
 
   int getGatingIndex() override { return startingIndexGating; };
 
+  void RHSUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
 };
 
 class MTensionModel : public MCellModel {
@@ -175,6 +173,10 @@ public:
   vector<double> YInfty(const vector<double> &var) const override {
     return elphy->YInfty(var);
   }
+
+  void RHSUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+
+  virtual void TensionUpdate(const vector<double> &var, vector<double> &g, double iExt) const {};
 };
 
 
diff --git a/src/cellmodels/MFitzHughNagumo.cpp b/src/cellmodels/MFitzHughNagumo.cpp
index 993ffdc61eb86e998e3e90d09395735f6c86738f..51d0db16dd5b87989a3e9c80550d09b3717cdae3 100644
--- a/src/cellmodels/MFitzHughNagumo.cpp
+++ b/src/cellmodels/MFitzHughNagumo.cpp
@@ -12,14 +12,6 @@ void MFitzHughNagumo::Initialize(vector<double> &var) {
   var[1] = 0.0;
 }
 
-void
-MFitzHughNagumo::EvaluateRHS(double t, const std::vector<double> &var,
-                             std::vector<double> &g) const {
-  double mVolt = var[0];
-  g[0] = c1 * mVolt * (mVolt - a) * (1 - mVolt) - c2 * var[1] + IExt(t);
-  g[1] = b * (mVolt - c3 * var[1]);
-}
-
 vector<double> MFitzHughNagumo::Tau(const vector<double> &var) const {
   return {};
 }
diff --git a/src/cellmodels/MFitzHughNagumo.hpp b/src/cellmodels/MFitzHughNagumo.hpp
index e6e932333a78d3d9327116755ffa214564b991c8..87276e5591757ef81f1c2328e6a30f9fee0d3a2b 100644
--- a/src/cellmodels/MFitzHughNagumo.hpp
+++ b/src/cellmodels/MFitzHughNagumo.hpp
@@ -28,8 +28,6 @@ public:
 
   void UpdateValues(double dt, vector<double> &var, const vector<double> &g) const override;
 
-  void EvaluateRHS(double t, const std::vector<double> &var, std::vector<double> &g) const override;
-
   double InitialValue() const override { return 0.0; };
 
   vector<double> Tau(const vector<double> &var) const override;
diff --git a/src/cellmodels/MRossi.cpp b/src/cellmodels/MRossi.cpp
index f23ece455dd6df2abaa83ddd36ffb6aed04f8d64..683ec1abe69f7d2223ae7726ab9d0b9e37f7f786 100644
--- a/src/cellmodels/MRossi.cpp
+++ b/src/cellmodels/MRossi.cpp
@@ -15,37 +15,6 @@ double MRossi::activeStressEvolution(double l, int N) const {
   return taylorSum;
 }
 
-void
-MRossi::EvaluateTensionRHS(double t, const std::vector<double> &var, std::vector<double> &g,
-                           int shift) const {
-  double gamma = var[shift];
-  double ca = var[elphy->CalciumIndex()] * 1e6;
-  double iota = var[shift + 2];
-
-  double activeForce =
-      -abs((ca > restingCA) *
-           sarcomereForce * (ca - restingCA) * (ca - restingCA) *
-           forceLengthRelation(sqrt(iota)));
-
-
-  double initialStress = iota * activeStressEvolution(gamma);
-
-  //mout << OUT(ca) << endl << activeForce + initialStress << endl;
-  double cubicLength = (1.0 + gamma) * (1.0 + gamma) * (1.0 + gamma);
-  double gshift = (0.001 / (viscosity * ca * ca)) *
-                  (activeForce + initialStress);
-
-/*
-  if(ca>restingCA){
-    mout << "iota = " << iota << " - sFa = " << sarcomereForce << " - RFL = " << forceLengthRelation(sqrt(iota)) << " - ";
-    mout << "ca = " << ca << " -  g = " << gshift << " - k = " << activeForce << endl;
-    mout << "Taylorexp = " << initialStress << endl;
-  }
-*/
-  g[shift] = gshift;// (2. *iota / cubicLength) - 2. * iota);
-  g[shift + 1] = activeForce;
-}
-
 double MRossi::forceLengthRelation(double stretch) const {
   //return (1-exp(-0.8*stretch));
 
@@ -64,8 +33,8 @@ double MRossi::forceLengthRelation(double stretch) const {
 
 void MRossi::Initialize(Vectors &vectors) {
   elphy->Initialize(vectors);
-  vectors[elphy->Size()] = 0.0;
-  vectors[elphy->Size() + 1] = 0.0;
+  R_gamma(vectors) = 0.0;
+  R_force(vectors) = 0.0;
 
   std::vector<double> elphyValues(vectors.size());
   elphy->Initialize(elphyValues);
@@ -74,13 +43,43 @@ void MRossi::Initialize(Vectors &vectors) {
 
 void MRossi::Initialize(vector<double> &vw) {
   elphy->Initialize(vw);
-  vw[elphy->Size()] = 0.0;
-  vw[elphy->Size() + 1] = 0.0;
+  R_gamma(vw) = 0.0;
+  R_force(vw) = 0.0;
 
   restingCA = vw[elphy->CalciumIndex()] * 1e6;
 }
 
-void MRossi::EvaluateRHS(double t, const vector<double> &var, vector<double> &g) const {
-  elphy->EvaluateRHS(t, var, g);
-  EvaluateTensionRHS(t, var, g, elphy->Size());
-}
\ No newline at end of file
+void MRossi::TensionUpdate(const vector<double> &var, vector<double> &g, double iExt) const {
+  double gamma = R_gamma(var);
+  double ca = var[elphy->CalciumIndex()] * 1e6;
+  double iota = R_iota(var);
+
+  double activeForce =
+      -abs((ca > restingCA) *
+           sarcomereForce * (ca - restingCA) * (ca - restingCA) *
+           forceLengthRelation(sqrt(iota)));
+
+
+  double initialStress = iota * activeStressEvolution(gamma);
+
+  //mout << OUT(ca) << endl << activeForce + initialStress << endl;
+  double cubicLength = (1.0 + gamma) * (1.0 + gamma) * (1.0 + gamma);
+  double gshift = (0.001 / (viscosity * ca * ca)) *
+                  (activeForce + initialStress);
+
+/*
+  if(ca>restingCA){
+    mout << "iota = " << iota << " - sFa = " << sarcomereForce << " - RFL = " << forceLengthRelation(sqrt(iota)) << " - ";
+    mout << "ca = " << ca << " -  g = " << gshift << " - k = " << activeForce << endl;
+    mout << "Taylorexp = " << initialStress << endl;
+  }
+*/
+  R_gamma(g) = gshift;// (2. *iota / cubicLength) - 2. * iota);
+  R_force(g) = activeForce;
+}
+
+void MRossi::UpdateValues(double dt, vector<double> &var, const vector<double> &g) const {
+  elphy->UpdateValues(dt, var, g);
+  R_gamma(var) += dt * R_gamma(g);
+  R_force(var) += dt * R_force(g);
+}
diff --git a/src/cellmodels/MRossi.hpp b/src/cellmodels/MRossi.hpp
index f12608f0f308147087de0b5a6f380d5104b8a2e9..e1693da5911ef85bcab784d82d8683799c9c3057 100644
--- a/src/cellmodels/MRossi.hpp
+++ b/src/cellmodels/MRossi.hpp
@@ -4,6 +4,11 @@
 
 #include "MCellModel.hpp"
 
+// Variable indexing
+#define R_gamma(v) v[elphy->Size()]
+#define R_force(v) v[elphy->Size()+1]
+#define R_iota(v) v[elphy->Size()+2]
+
 class MRossi : public MTensionModel {
 
   double restingCA{0.0};        // mM
@@ -24,36 +29,38 @@ class MRossi : public MTensionModel {
 
   double forceLengthRelation(double iota4) const;
 
+  double activeStressEvolution(double l, int N = 5) const;
+
 public:
-  explicit MRossi(std::unique_ptr<MElphyModel> &&M, double initCA=0.0) : MTensionModel(std::move(M), 1),
-                                                                         restingCA(initCA) {
+  explicit MRossi(std::unique_ptr<MElphyModel> &&M, double initCA = 0.0) : MTensionModel(
+      std::move(M), 1),
+                                                                           restingCA(initCA) {
     Config::Get("RossiViscosity", viscosity);
     Config::Get("RossiForce", sarcomereForce);
   }
 
+  double InitialValue() const override { return 0.0; }
+
   void Initialize(Vectors &VW) override;
 
   void Initialize(std::vector<double> &vw) override;
+
   /*
    * The vector var contains all intermediate variables in the following order:
    * 0. stretch
    * 1. calcium
    * 2. I4 (macroscopic second material invariant)
    */
-  void EvaluateTensionRHS(double t, const vector<double> &var, vector<double> &g,
-                          int shift) const override;
+  void TensionUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
 
-  double InitialValue() const override { return 0.0; }
+  void UpdateValues(double dt, vector<double> &var, const vector<double> &g) const override;
 
   int TensionIndex() const override { return elphy->Size() + 1; }
 
   int Size() const override { return elphy->Size() + 2; }
 
-  void EvaluateRHS(double t, const vector<double> &var, vector<double> &g) const override;
-
   double heavyside(double calcium);
 
-  double activeStressEvolution(double l, int N=5) const;
 };
 
 #endif //MROSSI_HPP
diff --git a/src/cellmodels/MTenTusscher.cpp b/src/cellmodels/MTenTusscher.cpp
index 9c6093d535609a35826399acce7de7d1adf87d4b..aba1469a3be6c45a44ad30c481e7a261e021189b 100644
--- a/src/cellmodels/MTenTusscher.cpp
+++ b/src/cellmodels/MTenTusscher.cpp
@@ -5,47 +5,47 @@ int MTenTusscher::CalciumIndex() const {
 }
 
 void MTenTusscher::Initialize(Vectors &VW) {
-  TT_Pot(VW) = V_init; //V // x
-  TT_Ca(VW) = 0.00007;//0.000153; //Ca // x
-  TT_Ca_SS(VW) = 0.00007;//0.00042; // Ca_SS // x
-  TT_Ca_SR(VW) = 1.3;//4.272; //Ca_SR // x
-  TT_Na(VW) = 7.67;//10.132; //Na // x
-  TT_K(VW) = 138.3;//138.52; //K // x
-  TT_m(VW) = 0.0; //0.00165; // m  // x
-  TT_h(VW) = 0.75; // 0.749; // h // x
-  TT_j(VW) = 0.75; // 0.6788; // j // x
-  TT_x_r1(VW) = 0.0; // 0.0165; // x_r1 // x
-  TT_x_r2(VW) = 1.0; // 0.473; // x_r2 // x
-  TT_x_s(VW) = 0.0; // 0.0174; // x_s // x
-  TT_s(VW) = 1.0; // 0.999998; // s // x
-  TT_r(VW) = 0.0; // 2.347e-8; // r // x
-  TT_d(VW) = 0.0; // 3.288e-5; // d // x
-  TT_f(VW) = 1.0; // 0.7026; // f // x
-  TT_f_2(VW) = 1.0; // 0.9526; // f_2 // x
-  TT_f_Ca(VW) = 1.0; // 0.9942; // f_Ca // x (falls f_cass)
-  TT_R_quer(VW) = 1.0; // 0.8978; // R_quer // x
+  TT_Pot(VW) = InitialValue(); //V // x
+  TT_Ca(VW) = 0.00007;//0.000153;
+  TT_Ca_SS(VW) = 0.00007;//0.00042;
+  TT_Ca_SR(VW) = 1.3;//4.272;
+  TT_Na(VW) = 7.67;//10.132;
+  TT_K(VW) = 138.3;//138.52;
+  TT_m(VW) = 0.0; //0.00165;
+  TT_h(VW) = 0.75; // 0.749;
+  TT_j(VW) = 0.75; // 0.6788;
+  TT_x_r1(VW) = 0.0; // 0.0165;
+  TT_x_r2(VW) = 1.0; // 0.473;
+  TT_x_s(VW) = 0.0; // 0.0174;
+  TT_s(VW) = 1.0; // 0.999998;
+  TT_r(VW) = 0.0; // 2.347e-8;
+  TT_d(VW) = 0.0; // 3.288e-5;
+  TT_f(VW) = 1.0; // 0.7026;
+  TT_f_2(VW) = 1.0; // 0.9526;
+  TT_f_Ca(VW) = 1.0; // 0.9942;
+  TT_R_quer(VW) = 1.0; // 0.8978;
 }
 
 void MTenTusscher::Initialize(vector<double> &vw) {
-  TT_Pot(vw) = V_init; //V // x
-  TT_Ca(vw) = 0.00007;//0.000153; //Ca // x
-  TT_Ca_SS(vw) = 0.00007;//0.00042; // Ca_SS // x
-  TT_Ca_SR(vw) = 1.3;//4.272; //Ca_SR // x
-  TT_Na(vw) = 7.67;//10.132; //Na // x
-  TT_K(vw) = 138.3;//138.52; //K // x
-  TT_m(vw) = 0.0;//0.00165;// m // x
-  TT_h(vw) = 0.75;//0.749; // h // x
-  TT_j(vw) = 0.75;//0.6788; // j // x
-  TT_x_r1(vw) = 0.0;//0.0165; // x_r1 // x
-  TT_x_r2(vw) = 1.0;//0.473; // x_r2 // x
-  TT_x_s(vw) = 0.0;//0.0174; // x_s // x
-  TT_s(vw) = 1.0;//0.999998; // s // x
-  TT_r(vw) = 0.0;//2.347e-8; // r // x
-  TT_d(vw) = 0.0;//3.288e-5; // d // x
-  TT_f(vw) = 1.0;//0.7026; // f // x
-  TT_f_2(vw) = 1.0;//0.9526; // f_2 // x
-  TT_f_Ca(vw) = 1.0;//0.9942; // f_Ca // x
-  TT_R_quer(vw) = 1.0;//0.8978; // R_quer // x
+  TT_Pot(vw) = InitialValue();
+  TT_Ca(vw) = 0.00007;//0.000153;
+  TT_Ca_SS(vw) = 0.00007;//0.00042;
+  TT_Ca_SR(vw) = 1.3;//4.272;
+  TT_Na(vw) = 7.67;//10.132;
+  TT_K(vw) = 138.3;//138.52;
+  TT_m(vw) = 0.0;//0.00165;
+  TT_h(vw) = 0.75;//0.749;
+  TT_j(vw) = 0.75;//0.6788;
+  TT_x_r1(vw) = 0.0;//0.0165;
+  TT_x_r2(vw) = 1.0;//0.473;
+  TT_x_s(vw) = 0.0;//0.0174;
+  TT_s(vw) = 1.0;//0.999998;
+  TT_r(vw) = 0.0;//2.347e-8;
+  TT_d(vw) = 0.0;//3.288e-5;
+  TT_f(vw) = 1.0;//0.7026;
+  TT_f_2(vw) = 1.0;//0.9526;
+  TT_f_Ca(vw) = 1.0;//0.9942;
+  TT_R_quer(vw) = 1.0;//0.8978;
 }
 
 double tauFormula_h(double V) {
diff --git a/src/cellmodels/MTenTusscher.hpp b/src/cellmodels/MTenTusscher.hpp
index 8b56f57b067a621c81da0b1992bbe44eff6414d1..6b2bec600bc6d3332e30c947a6102f8d6fc2cfa8 100644
--- a/src/cellmodels/MTenTusscher.hpp
+++ b/src/cellmodels/MTenTusscher.hpp
@@ -84,7 +84,6 @@ class MTenTusscher : public MElphyModel {
   double max_sr = 2.5;
   double min_sr = 1.0;
   double EC = 1.5;
-  double V_init= -86.2;//-85.423;
 
   //Calcium buffering dynamics
   double Bufc=0.2;
@@ -112,19 +111,22 @@ public:
   int CalciumIndex() const override;
 
 
+  double InitialValue() const override { return -86.2; } //-85.423
+
   void Initialize(Vectors &VW) override;
 
   void Initialize(std::vector<double> &vw) override;
 
 
   void IIonUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
-  void ConcentrationUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+
+  void
+  ConcentrationUpdate(const vector<double> &var, vector<double> &g, double iExt) const override;
+
   void GatingUpdate(const vector<double> &var, vector<double> &g) const override;
-  void ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const override;
 
+  void ExpIntegratorUpdate(double dt, const vector<double> &var, vector<double> &g) const override;
 
-  void EvaluateRHS(double t, const std::vector<double> &var, std::vector<double> &g) const override{};
-  double InitialValue() const override { return V_init; };
 
   vector<double> Tau(const vector<double> &var) const override;
   vector<double> YInfty(const vector<double> &var) const override;
diff --git a/src/cellmodels/solvers/MElphySolver.cpp b/src/cellmodels/solvers/MElphySolver.cpp
index c77aa5ba146da420c6296e0be54a4c167694bdc1..3951d4b487debccdfa4b0b040f075e151fd90097 100644
--- a/src/cellmodels/solvers/MElphySolver.cpp
+++ b/src/cellmodels/solvers/MElphySolver.cpp
@@ -249,7 +249,7 @@ void EIWithNewtonStep(MCellModel &cellModel, double t, double dt, vector<double>
   int position = cellModel.getGatingIndex();
 
   cellModel.ExpIntegratorUpdate(dt, VW, g[0]);
-  cellModel.EvaluateRHS(t, VW, g[0]);
+  cellModel.RHSUpdate(VW, g[0], cellModel.IExt(t));
   double DVIion = cellModel.JacobiEntry(0, 0, VW);
   VW[0] = (1 / (1 + dt * DVIion)) * ((1 + dt * DVIion) * mVolt + dt * g[0][0]);
   cellModel.ConcentrationUpdate(VW, g[0], cellModel.IExt(t));
@@ -269,7 +269,7 @@ void EigenvalueStep(MCellModel &cellModel, double t, double dt, vector<double> &
     double tau_yInverse = alphaV + betaV;
     VW[i] = y_infty + (VW[i] - y_infty) * exp(-1000 * dt * tau_yInverse);
   }
-  cellModel.EvaluateRHS(t, VW, g[0]);
+  cellModel.RHSUpdate(VW, g[0], cellModel.IExt(t));
   VW[0] += dt * 1000.0 * g[0][0];
   VW[1] += dt * 1000.0 * g[0][1];
 }
diff --git a/test/cellmodels/MTest.hpp b/test/cellmodels/MTest.hpp
index 9db83831190acc5e81e521ebc2524d23b0357737..6a3b9ac61cd81da0cd9fc79455fb059728cea91b 100644
--- a/test/cellmodels/MTest.hpp
+++ b/test/cellmodels/MTest.hpp
@@ -8,8 +8,6 @@ public:
   explicit MTest() : MElphyModel(0, 0,1.0) {
   }
 
-  CellModelType Type() const override { return ELPHY; }
-
   void Initialize(Vectors &VW) override {
     VW[0] = 0;
   };
@@ -18,11 +16,18 @@ public:
     InitialValues(0, vw);
   }
 
-  void
-  EvaluateRHS(double t, const std::vector<double> &var, std::vector<double> &g) const override {
+  double IExt(double t) const override {
+    return t;
+  }
+
+  void RHSUpdate(const vector<double> &var, vector<double> &g, double t) const override {
     g[0] = Diff(t, var);
   }
 
+  void UpdateValues(double dt, vector<double> &var, const vector<double> &g) const override {
+    var[0] += dt * g[0];
+  }
+
   virtual void InitialValues(double t, std::vector<double> &vw) = 0;
 
   virtual double Diff(double t, const std::vector<double> &vw) const = 0;