diff --git a/3rdParty b/3rdParty
index 0b3ca6f45..fbcbfc76f 160000
--- a/3rdParty
+++ b/3rdParty
@@ -1 +1 @@
-Subproject commit 0b3ca6f4576c6b1bde14bc44658960fc1d78de15
+Subproject commit fbcbfc76f3f71617ba6f006c0603c0049ebce5a2
diff --git a/src/OMSimulatorLib/Component.h b/src/OMSimulatorLib/Component.h
index cac09a62c..122bee289 100644
--- a/src/OMSimulatorLib/Component.h
+++ b/src/OMSimulatorLib/Component.h
@@ -46,6 +46,7 @@
 #include "ResultWriter.h"
 #include "SignalDerivative.h"
 #include "Snapshot.h"
+#include "Variable.h"
 #include "OMSimulator/Types.h"
 #include "Values.h"
 
@@ -79,6 +80,14 @@ namespace oms
     virtual oms_status_enu_t terminate() = 0;
     virtual oms_status_enu_t updateSignals(ResultWriter& resultWriter) = 0;
     virtual Variable* getVariable(const ComRef& cref) = 0;
+    virtual oms_status_enu_t getVariableType(const ComRef& cref, oms_signal_type_enu_t& type)
+    {
+      Variable* var = getVariable(cref);
+      if (!var)
+        return oms_status_error;
+      type = var->getType();
+      return oms_status_ok;
+    }
 
     virtual bool getCanGetAndSetState() { return false; }
     virtual const FMUInfo* getFMUInfo() const { return nullptr; }
diff --git a/src/OMSimulatorLib/ComponentFMU3CS.cpp b/src/OMSimulatorLib/ComponentFMU3CS.cpp
index b6e41cdd8..3ed39de58 100644
--- a/src/OMSimulatorLib/ComponentFMU3CS.cpp
+++ b/src/OMSimulatorLib/ComponentFMU3CS.cpp
@@ -688,7 +688,7 @@ oms_status_enu_t oms::ComponentFMU3CS::instantiate()
   time = getModel().getStartTime();
 
   double relativeTolerance = 0.0;
-  dynamic_cast<SystemWC*>(getParentSystem())->getTolerance(&relativeTolerance);
+  getParentSystem()->getTolerance(&relativeTolerance);
 
   fmi3Status status_ = fmi3_enterInitializationMode(fmu, fmi3False, relativeTolerance, time, fmi3False, getModel().getStopTime());
 
@@ -882,7 +882,7 @@ oms_status_enu_t oms::ComponentFMU3CS::reset()
   // enterInitialization
   time = getModel().getStartTime();
   double relativeTolerance = 0.0;
-  dynamic_cast<SystemWC*>(getParentSystem())->getTolerance(&relativeTolerance);
+  getParentSystem()->getTolerance(&relativeTolerance);
 
   fmi3Status status_ = fmi3_enterInitializationMode(fmu, fmi3False, relativeTolerance, time, fmi3True, getModel().getStopTime());
   if (fmi3OK != fmistatus) return logError_FMUCall("fmi3_enterInitializationMode", this);
diff --git a/src/OMSimulatorLib/ComponentFMU3ME.cpp b/src/OMSimulatorLib/ComponentFMU3ME.cpp
index d07ee7725..31b1ff644 100644
--- a/src/OMSimulatorLib/ComponentFMU3ME.cpp
+++ b/src/OMSimulatorLib/ComponentFMU3ME.cpp
@@ -675,7 +675,7 @@ oms_status_enu_t oms::ComponentFMU3ME::instantiate()
   const double& startTime = getModel().getStartTime();
 
   double relativeTolerance = 0.0;
-  dynamic_cast<SystemSC*>(getParentSystem())->getTolerance(&relativeTolerance);
+  getParentSystem()->getTolerance(&relativeTolerance);
 
   fmi3Status status_ = fmi3_enterInitializationMode(fmu, fmi3False, relativeTolerance, startTime, fmi3False, getModel().getStopTime());
 
@@ -913,7 +913,7 @@ oms_status_enu_t oms::ComponentFMU3ME::reset()
   // enterInitialization
   const double& startTime = getModel().getStartTime();
   double relativeTolerance = 0.0;
-  dynamic_cast<SystemSC3*>(getParentSystem())->getTolerance(&relativeTolerance);
+  getParentSystem()->getTolerance(&relativeTolerance);
 
   fmi3Status status_ = fmi3_enterInitializationMode(fmu, fmi3False, relativeTolerance, startTime, fmi3True, getModel().getStopTime());
   if (fmi3OK != fmistatus) return logError_FMUCall("fmi3_enterInitializationMode", this);
diff --git a/src/OMSimulatorLib/ComponentFMUCS.cpp b/src/OMSimulatorLib/ComponentFMUCS.cpp
index 3b301bee1..2b9715de4 100644
--- a/src/OMSimulatorLib/ComponentFMUCS.cpp
+++ b/src/OMSimulatorLib/ComponentFMUCS.cpp
@@ -674,7 +674,7 @@ oms_status_enu_t oms::ComponentFMUCS::instantiate()
   // enterInitialization
   time = getModel().getStartTime();
   double relativeTolerance = 0.0;
-  dynamic_cast<SystemWC*>(getParentSystem())->getTolerance(&relativeTolerance);
+  getParentSystem()->getTolerance(&relativeTolerance);
 
   fmi2Status status = fmi2_setupExperiment(fmu, fmi2True, relativeTolerance, time, fmi2False, 1.0);
   if (fmi2OK != status) return logError_FMUCall("fmi2_setupExperiment", this);
@@ -868,7 +868,7 @@ oms_status_enu_t oms::ComponentFMUCS::reset()
   // enterInitialization
   time = getModel().getStartTime();
   double relativeTolerance = 0.0;
-  dynamic_cast<SystemWC*>(getParentSystem())->getTolerance(&relativeTolerance);
+  getParentSystem()->getTolerance(&relativeTolerance);
   fmistatus = fmi2_setupExperiment(fmu, fmi2True, relativeTolerance, time, fmi2False, 1.0);
   if (fmi2OK != fmistatus) return logError_FMUCall("fmi2_setupExperiment", this);
 
diff --git a/src/OMSimulatorLib/ComponentFMUME.cpp b/src/OMSimulatorLib/ComponentFMUME.cpp
index eb929a153..b31042fcc 100644
--- a/src/OMSimulatorLib/ComponentFMUME.cpp
+++ b/src/OMSimulatorLib/ComponentFMUME.cpp
@@ -48,7 +48,7 @@
 #include <regex>
 #include <unordered_set>
 #include <cmath>
-
+#include <iostream>
 oms::ComponentFMUME::ComponentFMUME(const ComRef& cref, System* parentSystem, const std::string& fmuPath)
   : oms::Component(cref, oms_component_fmu, parentSystem, fmuPath), fmuInfo(fmuPath)
 {
@@ -704,7 +704,8 @@ oms_status_enu_t oms::ComponentFMUME::instantiate()
   // enterInitialization
   const double& startTime = getModel().getStartTime();
   double relativeTolerance = 0.0;
-  dynamic_cast<SystemSC*>(getParentSystem())->getTolerance(&relativeTolerance);
+  getParentSystem()->getTolerance(&relativeTolerance);
+
   fmi2Status status = fmi2_setupExperiment(fmu, fmi2True, relativeTolerance, startTime, fmi2False, 1.0);
   if (fmi2OK != status) return logError_FMUCall("fmi2_setupExperiment", this);
 
@@ -929,7 +930,7 @@ oms_status_enu_t oms::ComponentFMUME::reset()
   // enterInitialization
   const double& startTime = getModel().getStartTime();
   double relativeTolerance = 0.0;
-  dynamic_cast<SystemSC*>(getParentSystem())->getTolerance(&relativeTolerance);
+  getParentSystem()->getTolerance(&relativeTolerance);
   fmistatus = fmi2_setupExperiment(fmu, fmi2True, relativeTolerance, startTime, fmi2False, 1.0);
   if (fmi2OK != fmistatus) return logError_FMUCall("fmi2_setupExperiment", this);
 
diff --git a/src/OMSimulatorLib/ComponentTable.h b/src/OMSimulatorLib/ComponentTable.h
index c4219ab53..1c5348526 100644
--- a/src/OMSimulatorLib/ComponentTable.h
+++ b/src/OMSimulatorLib/ComponentTable.h
@@ -70,6 +70,8 @@ namespace oms
     oms_status_enu_t setTime(double time) {this->time = time; return oms_status_ok;}
 
     Variable* getVariable(const ComRef& cref) {logError_NotImplemented; return NULL;}
+    // All ComponentTable signals are interpolated series values (see getReal/getInteger/getBoolean below).
+    oms_status_enu_t getVariableType(const ComRef& cref, oms_signal_type_enu_t& type) {type = oms_signal_type_real; return oms_status_ok;}
     oms_status_enu_t getReal(const ComRef& cref, double& value);
     oms_status_enu_t getRealOutputDerivative(const ComRef& cref, SignalDerivative& value);
     oms_status_enu_t getInteger(const ComRef& cref, int& value);
diff --git a/src/OMSimulatorLib/System.cpp b/src/OMSimulatorLib/System.cpp
index a68f7b770..32c783577 100644
--- a/src/OMSimulatorLib/System.cpp
+++ b/src/OMSimulatorLib/System.cpp
@@ -204,11 +204,7 @@ oms_status_enu_t oms::System::getVariableType(const ComRef& cref, oms_signal_typ
 
   auto component = components.find(head);
   if (component != components.end())
-  {
-    oms::Variable * var = component->second->getVariable(tail);
-    value = var->getType();
-    return oms_status_ok;
-  }
+    return component->second->getVariableType(tail, value);
 
   // check connectors
   for (auto& connector : connectors)
diff --git a/src/OMSimulatorLib/SystemSC.cpp b/src/OMSimulatorLib/SystemSC.cpp
index c7cfc5d0d..1818d4c8c 100644
--- a/src/OMSimulatorLib/SystemSC.cpp
+++ b/src/OMSimulatorLib/SystemSC.cpp
@@ -52,10 +52,9 @@ int oms::cvode_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data)
 {
   SystemSC* system = (SystemSC*)user_data;
   oms_status_enu_t status;
-  fmi2Status fmistatus;
 
-  // update states in FMUs
-  for (size_t i=0, j_y=0, j_ydot=0; i < system->fmus.size(); ++i)
+  // set time and states in all FMUs for this residual evaluation
+  for (size_t i=0, j_y=0; i < system->fmus.size(); ++i)
   {
     system->fmus[i]->setTime(t);
 
@@ -68,8 +67,25 @@ int oms::cvode_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data)
     // set states
     status = system->fmus[i]->setContinuousStates(system->states[i]);
     if (oms_status_ok != status) return status;
+  }
+
+  // Propagate connected outputs into inputs (e.g. algebraic feedthrough across
+  // FMU boundaries) so that the derivatives below see coupling inputs that are
+  // consistent with the current state guess, instead of the value left over
+  // from the last accepted step. Without this, tightly-coupled FMUs (e.g. two
+  // FMUs closing an algebraic loop through connected signals) only get their
+  // coupling updated once per accepted step, which is an explicit one-step-
+  // delayed (Jacobi-style) coupling hidden inside what should be a fully
+  // implicit combined ODE solve, and can destabilize the integration.
+  status = system->updateInputs(system->simulationGraph);
+  if (oms_status_ok != status) return status;
+
+  // get state derivatives
+  for (size_t i=0, j_ydot=0; i < system->fmus.size(); ++i)
+  {
+    if (0 == system->nStates[i])
+      continue;
 
-    // get state derivatives
     status = system->fmus[i]->getDerivatives(system->states_der[i]);
     if (oms_status_ok != status) return status;
 
@@ -537,6 +553,25 @@ oms_status_enu_t oms::SystemSC::reset()
     solverData.cvode.mem = nullptr;
   }
 
+  // Free arrays allocated by initialize() and clear the per-FMU tracking
+  // vectors it populated, so a subsequent initialize() call starts fresh
+  // instead of appending duplicate entries for every FMU (which corrupts
+  // solver setup and causes lifecycle calls like enterContinuousTimeMode()
+  // to be issued twice for the same FMU).
+  for (double* ptr : states) free(ptr);
+  for (double* ptr : states_der) free(ptr);
+  for (double* ptr : states_nominal) free(ptr);
+  for (double* ptr : event_indicators) free(ptr);
+  for (double* ptr : event_indicators_prev) free(ptr);
+  fmus.clear();
+  nStates.clear();
+  nEventIndicators.clear();
+  states.clear();
+  states_der.clear();
+  states_nominal.clear();
+  event_indicators.clear();
+  event_indicators_prev.clear();
+
   return oms_status_ok;
 }
 
diff --git a/src/OMSimulatorLib/SystemSC3.cpp b/src/OMSimulatorLib/SystemSC3.cpp
index c4cd4862e..87c9fdf0c 100644
--- a/src/OMSimulatorLib/SystemSC3.cpp
+++ b/src/OMSimulatorLib/SystemSC3.cpp
@@ -51,10 +51,9 @@ int oms::cvode_rhs3(realtype t, N_Vector y, N_Vector ydot, void* user_data)
 {
   SystemSC3* system = (SystemSC3*)user_data;
   oms_status_enu_t status;
-  fmi3Status fmistatus;
 
-  // update states in FMUs
-  for (size_t i=0, j_y=0, j_ydot=0; i < system->fmus.size(); ++i)
+  // set time and states in all FMUs for this residual evaluation
+  for (size_t i=0, j_y=0; i < system->fmus.size(); ++i)
   {
     system->fmus[i]->setTime(t);
 
@@ -67,8 +66,24 @@ int oms::cvode_rhs3(realtype t, N_Vector y, N_Vector ydot, void* user_data)
     // set states
     status = system->fmus[i]->setContinuousStates(system->states[i]);
     if (oms_status_ok != status) return status;
+  }
+
+  // Propagate connected outputs into inputs (e.g. algebraic feedthrough across
+  // FMU boundaries) so that the derivatives below see coupling inputs that are
+  // consistent with the current state guess, instead of the value left over
+  // from the last accepted step. Without this, tightly-coupled FMUs only get
+  // their coupling updated once per accepted step, which is an explicit
+  // one-step-delayed (Jacobi-style) coupling hidden inside what should be a
+  // fully implicit combined ODE solve, and can destabilize the integration.
+  status = system->updateInputs(system->simulationGraph);
+  if (oms_status_ok != status) return status;
+
+  // get state derivatives
+  for (size_t i=0, j_ydot=0; i < system->fmus.size(); ++i)
+  {
+    if (0 == system->nStates[i])
+      continue;
 
-    // get state derivatives
     status = system->fmus[i]->getDerivatives(system->states_der[i]);
     if (oms_status_ok != status) return status;
 
@@ -531,6 +546,25 @@ oms_status_enu_t oms::SystemSC3::reset()
     solverData.cvode.mem = nullptr;
   }
 
+  // Free arrays allocated by initialize() and clear the per-FMU tracking
+  // vectors it populated, so a subsequent initialize() call starts fresh
+  // instead of appending duplicate entries for every FMU (which corrupts
+  // solver setup and causes lifecycle calls like enterContinuousTimeMode()
+  // to be issued twice for the same FMU).
+  for (double* ptr : states) free(ptr);
+  for (double* ptr : states_der) free(ptr);
+  for (double* ptr : states_nominal) free(ptr);
+  for (double* ptr : event_indicators) free(ptr);
+  for (double* ptr : event_indicators_prev) free(ptr);
+  fmus.clear();
+  nStates.clear();
+  nEventIndicators.clear();
+  states.clear();
+  states_der.clear();
+  states_nominal.clear();
+  event_indicators.clear();
+  event_indicators_prev.clear();
+
   return oms_status_ok;
 }
 
diff --git a/src/OMSimulatorPython/OMSimulatorPython3.bat b/src/OMSimulatorPython/OMSimulatorPython3.bat
index 5ebeb46ab..e9257fba9 100644
--- a/src/OMSimulatorPython/OMSimulatorPython3.bat
+++ b/src/OMSimulatorPython/OMSimulatorPython3.bat
@@ -1,6 +1,6 @@
 @echo off
 
-set PYTHONPATH=%~dp0;%~dp0\..\lib;%PYTHONPATH%
+set PYTHONPATH=%~dp0;%~dp0\..\lib;%~dp0\..\lib\omc;%PYTHONPATH%
 set PATH=%~dp0;%~dp0\..\bin;%PATH%
 
 python3 %*
diff --git a/src/OMSimulatorPython/OMSimulatorPython3.in b/src/OMSimulatorPython/OMSimulatorPython3.in
index 11c2c22ec..9af5a2663 100644
--- a/src/OMSimulatorPython/OMSimulatorPython3.in
+++ b/src/OMSimulatorPython/OMSimulatorPython3.in
@@ -6,7 +6,7 @@ SCRIPTPATH=$(dirname "$SCRIPT")
 OMS_LIB_DIR="${SCRIPTPATH}/../@CMAKE_INSTALL_LIBDIR@"
 
 # Add the oms install lib directory to PYTHONPATH.
-if [ -z "$PYTHONPATH" ]; then export PYTHONPATH="${OMS_LIB_DIR}"; else export PYTHONPATH="${OMS_LIB_DIR}":$PYTHONPATH; fi
+if [ -z "$PYTHONPATH" ]; then export PYTHONPATH="${OMS_LIB_DIR}:${OMS_LIB_DIR}/omc"; else export PYTHONPATH="${OMS_LIB_DIR}:${OMS_LIB_DIR}/omc":$PYTHONPATH; fi
 
 # On Linux, add the oms install LIB dir to LD_LIBRARY_PATH. That is where the shard lib is installed to.
 # On Windows, add the oms install BIN dir to PATH. That is where the dll is installed to.
diff --git a/src/OMSimulatorPython/capi.py b/src/OMSimulatorPython/capi.py
index 4be89da00..38b8a88f0 100644
--- a/src/OMSimulatorPython/capi.py
+++ b/src/OMSimulatorPython/capi.py
@@ -106,6 +106,8 @@ def __init__(self):
     self.obj.oms_newModel.restype = ctypes.c_int
     self.obj.oms_removeSignalsFromResults.argtypes = [ctypes.c_char_p, ctypes.c_char_p]
     self.obj.oms_removeSignalsFromResults.restype = ctypes.c_int
+    self.obj.oms_reset.argtypes = [ctypes.c_char_p]
+    self.obj.oms_reset.restype = ctypes.c_int
     self.obj.oms_setCommandLineOption.argtypes = [ctypes.c_char_p]
     self.obj.oms_setCommandLineOption.restype = ctypes.c_int
     self.obj.oms_setConnectionLinearTransformation.argtypes = [ctypes.c_char_p, ctypes.c_char_p, ctypes.c_double, ctypes.c_double]
@@ -379,4 +381,9 @@ def removeSignalsFromResults(self, cref, regex):
     status = self.obj.oms_removeSignalsFromResults(cref.encode(), regex.encode())
     return Status(status)
 
+  def reset(self, cref) -> Status:
+    '''Reset the model to the state right after instantiation, discarding initialization/simulation progress.'''
+    status = self.obj.oms_reset(cref.encode())
+    return Status(status)
+
 Capi = capi()
diff --git a/src/OMSimulatorPython/componenttable.py b/src/OMSimulatorPython/componenttable.py
index 210c67df5..bc27d4117 100644
--- a/src/OMSimulatorPython/componenttable.py
+++ b/src/OMSimulatorPython/componenttable.py
@@ -62,17 +62,12 @@ def __init__(self, name: CRef, filePath: Path | str, connectors : list | None =
     self.name = name
     self.filePath = filePath
     self.connectors = connectors
-    self._solver = None
+    self.solver = None
     self.elementgeometry = None
 
 
-  @property
-  def solver(self):
-    return self._solver
-
-  @solver.setter
-  def solver(self, cref: str):
-    self._solver = cref
+  def setSolver(self, cref: str):
+    self.solver = cref
 
   def list(self, prefix=""):
     print(f"{prefix} Table: {self.name}")
diff --git a/src/OMSimulatorPython/instantiated_model.py b/src/OMSimulatorPython/instantiated_model.py
index 8b3885de4..dbf830515 100644
--- a/src/OMSimulatorPython/instantiated_model.py
+++ b/src/OMSimulatorPython/instantiated_model.py
@@ -567,6 +567,12 @@ def initialize(self):
     if status != Status.ok:
       raise RuntimeError(f"Failed to initialize model: {status}")
 
+  def reset(self):
+    '''Reset the model to the state right after instantiation, discarding initialization/simulation progress.'''
+    status = Capi.reset(self.modelName)
+    if status != Status.ok:
+      raise RuntimeError(f"Failed to reset model: {status}")
+
   def simulate(self):
     status = Capi.simulate(self.modelName)
     if status != Status.ok:
diff --git a/testsuite/dcmotor/dcmotor.py b/testsuite/dcmotor/dcmotor.py
index e13b5ad9e..c913c9e08 100644
--- a/testsuite/dcmotor/dcmotor.py
+++ b/testsuite/dcmotor/dcmotor.py
@@ -212,8 +212,8 @@
 ## info:    maximum step size for 'model.root': 0.001000
 ## info:    Result file: dc_motor_res.mat (bufferSize=10)
 ## info:    Final Statistics for 'model.root':
-##          NumSteps = 10663 NumRhsEvals  = 10833 NumLinSolvSetups = 606
-##          NumNonlinSolvIters = 10832 NumNonlinSolvConvFails = 0 NumErrTestFails = 26
+##          NumSteps = 10030 NumRhsEvals  = 10044 NumLinSolvSetups = 511
+##          NumNonlinSolvIters = 10043 NumNonlinSolvConvFails = 0 NumErrTestFails = 1
 ## info:    3 warnings
 ## info:    0 errors
 ## endResult
diff --git a/testsuite/simulation/DualMassOscillator.py b/testsuite/simulation/DualMassOscillator1.py
similarity index 87%
rename from testsuite/simulation/DualMassOscillator.py
rename to testsuite/simulation/DualMassOscillator1.py
index d000df64d..fe716f620 100644
--- a/testsuite/simulation/DualMassOscillator.py
+++ b/testsuite/simulation/DualMassOscillator1.py
@@ -1,5 +1,5 @@
 ## status: correct
-## teardown_command: rm -rf DualMassOscillator.ssp DualMassOscillator_res.mat
+## teardown_command: rm -rf DualMassOscillator1.ssp DualMassOscillator1_res.mat
 ## linux: yes
 ## ucrt64: yes
 ## win: yes
@@ -24,12 +24,12 @@
 model.addConnection(CRef('default', 'system1', 'v1'), CRef('default', 'system2', 'v1'))
 model.addConnection(CRef('default', 'system1', 'a1'), CRef('default', 'system2', 'a1'))
 
-model.export('DualMassOscillator.ssp')
+model.export('DualMassOscillator1.ssp')
 
-model2 = SSP('DualMassOscillator.ssp')
+model2 = SSP('DualMassOscillator1.ssp')
 
 instantiated_model = model2.instantiate() ## internally generate the json file and also set the model state like virgin,
-instantiated_model.setResultFile("DualMassOscillator_res.mat")
+instantiated_model.setResultFile("DualMassOscillator1_res.mat")
 instantiated_model.setStopTime(10)
 
 instantiated_model.setFixedStepSize(1e-4)
@@ -52,7 +52,7 @@
 instantiated_model.delete()
 
 ## Result:
-## info:    Result file: DualMassOscillator_res.mat (bufferSize=10)
+## info:    Result file: DualMassOscillator1_res.mat (bufferSize=10)
 ## info:    Initialization
 ## info:      default.system1.x1: 0.0
 ## info:      default.system2.x2: 0.5
diff --git a/testsuite/simulation/DualMassOscillator2.py b/testsuite/simulation/DualMassOscillator2.py
new file mode 100644
index 000000000..124b7dcf9
--- /dev/null
+++ b/testsuite/simulation/DualMassOscillator2.py
@@ -0,0 +1,73 @@
+## status: correct
+## teardown_command: rm -rf DualMassOscillator2.ssp DualMassOscillator2_res.mat
+## linux: yes
+## ucrt64: yes
+## win: yes
+## mac: yes
+
+from OMSimulator import SSP, CRef, Settings
+
+Settings.suppressPath = True
+
+
+model = SSP()
+model.addResource('../resources/DualMassOscillator.System1.fmu', new_name='resources/system1.fmu')
+model.addResource('../resources/DualMassOscillator.System2.fmu', new_name='resources/system2.fmu')
+
+
+model.addComponent(CRef('default', 'system1'), 'resources/system1.fmu')
+model.addComponent(CRef('default', 'system2'), 'resources/system2.fmu')
+
+## add connection
+model.addConnection(CRef('default', 'system1', 'F'),  CRef('default', 'system2', 'F'))
+model.addConnection(CRef('default', 'system1', 'x1'), CRef('default', 'system2', 'x1'))
+model.addConnection(CRef('default', 'system1', 'v1'), CRef('default', 'system2', 'v1'))
+model.addConnection(CRef('default', 'system1', 'a1'), CRef('default', 'system2', 'a1'))
+
+solver1 = {'name' : 'solver1',  'method': 'cvode', 'relativeTolerance': 1e-6, 'initialStepSize': 1e-12, 'minimumStepSize': 1e-12, 'maximumStepSize': 1e-3}
+model.newSolver(solver1)
+
+model.setSolver(CRef('default', 'system1'), 'solver1')
+model.setSolver(CRef('default', 'system2'), 'solver1')
+
+
+model.export('DualMassOscillator2.ssp')
+
+model2 = SSP('DualMassOscillator2.ssp')
+
+instantiated_model = model2.instantiate() ## internally generate the json file and also set the model state like virgin,
+instantiated_model.setResultFile("DualMassOscillator2_res.mat")
+instantiated_model.setStopTime(10)
+
+instantiated_model.setFixedStepSize(1e-4)
+
+## set value
+instantiated_model.setValue(CRef('default', 'system1', 'x1_start'), 0.0)
+instantiated_model.setValue(CRef('default', 'system2', 'x2_start'), 0.5)
+
+instantiated_model.initialize()
+print(f"info:    Initialization")
+print(f"info:      default.system1.x1: {instantiated_model.getValue(CRef('default', 'system1', 'x1'))}")
+print(f"info:      default.system2.x2: {instantiated_model.getValue(CRef('default', 'system2', 'x2'))}")
+
+instantiated_model.simulate()
+print(f"info:    Simulation")
+print(f"info:      default.system1.x1: {instantiated_model.getValue(CRef('default', 'system1', 'x1'))}")
+print(f"info:      default.system2.x2: {instantiated_model.getValue(CRef('default', 'system2', 'x2'))}")
+
+instantiated_model.terminate()
+instantiated_model.delete()
+
+## Result:
+## info:    maximum step size for 'model.root': 0.000100
+## info:    Result file: DualMassOscillator2_res.mat (bufferSize=10)
+## info:    Final Statistics for 'model.root':
+##          NumSteps = 100001 NumRhsEvals  = 100002 NumLinSolvSetups = 5001
+##          NumNonlinSolvIters = 100001 NumNonlinSolvConvFails = 0 NumErrTestFails = 0
+## info:    Initialization
+## info:      default.system1.x1: 0.0
+## info:      default.system2.x2: 0.5
+## info:    Simulation
+## info:      default.system1.x1: 0.05067159474969803
+## info:      default.system2.x2: 0.03142555558313003
+## endResult
diff --git a/testsuite/simulation/Makefile b/testsuite/simulation/Makefile
index 5c1b71d38..09721c7cc 100644
--- a/testsuite/simulation/Makefile
+++ b/testsuite/simulation/Makefile
@@ -1,7 +1,8 @@
 TEST = ../rtest -v
 
 TESTFILES = \
-DualMassOscillator.py \
+DualMassOscillator1.py \
+DualMassOscillator2.py \
 exportJson1.py \
 exportJson2.py \
 exportJson4.py \