diff --git a/CMakeLists.txt b/CMakeLists.txt
index f6327bc..72ec21f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -47,6 +47,7 @@ if ( BarElm_LIBRARY )
 endif ( BarElm_LIBRARY )
 
 if ( BUILD_SOLVER_AS_DLL )
+  # Check if python is installed
   if ( CMAKE_VERSION VERSION_LESS 3.12 )
     find_package ( PythonInterp 3 ) # sets ${PYTHON_EXECUTABLE}
   else ( CMAKE_VERSION VERSION_LESS 3.12 )
@@ -60,6 +61,8 @@ if ( BUILD_SOLVER_AS_DLL )
   endif ( PythonInterp_FOUND )
 endif ( BUILD_SOLVER_AS_DLL )
 
+################################################################################
+
 set ( FT_TOLERANCE 1.0e-12 CACHE STRING "Regression test comparison tolerance" )
 option ( FT_LARGE_MODELS "Include large regression test models" false )
 option ( FT_RUN_VALGRIND "Run regression tests with valgrind" false )
@@ -77,15 +80,20 @@ if ( VALGRIND )
   set ( VALGRIND_CMD ${VALGRIND} --leak-check=full --show-leak-kinds=all )
 endif ( VALGRIND )
 
+################################################################################
+# Create a custom target to copy relevant input files into the build tree.
+
 add_custom_target ( copy_test_input_files )
 foreach ( TEST ${TEST_DIRS} )
+  set ( TIME_DOMAIN TimeDomain/${TEST} )
+  set ( FREQ_DOMAIN FrequencyDomain/${TEST} )
   set ( BIN_DIR ${CMAKE_CURRENT_BINARY_DIR}/${TEST} )
-  set ( CML_TDA ${CMAKE_CURRENT_SOURCE_DIR}/TimeDomain/${TEST}/CMakeLists.txt )
-  set ( CML_FDA ${CMAKE_CURRENT_SOURCE_DIR}/FrequencyDomain/${TEST}/CMakeLists.txt )
+  set ( CML_TDA ${CMAKE_CURRENT_SOURCE_DIR}/${TIME_DOMAIN}/CMakeLists.txt )
+  set ( CML_FDA ${CMAKE_CURRENT_SOURCE_DIR}/${FREQ_DOMAIN}/CMakeLists.txt )
   if ( EXISTS ${CML_TDA} )
-    set ( SRC_DIR TimeDomain/${TEST} )
+    set ( SRC_DIR ${TIME_DOMAIN} )
   elseif ( EXISTS ${CML_FDA} AND FFT_LIBRARY )
-    set ( SRC_DIR FrequencyDomain/${TEST} )
+    set ( SRC_DIR ${FREQ_DOMAIN} )
   elseif ( EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${TEST}/CMakeLists.txt )
     set ( SRC_DIR ${TEST} )
   else ( EXISTS ${CML_TDA} )
@@ -94,9 +102,17 @@ foreach ( TEST ${TEST_DIRS} )
   add_subdirectory ( ${SRC_DIR} ${BIN_DIR} )
   set ( SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/${SRC_DIR} )
   add_custom_command ( TARGET copy_test_input_files POST_BUILD
-                       COMMAND ${CMAKE_COMMAND} -E copy_directory ${SRC_DIR} ${BIN_DIR} )
+                       COMMAND ${CMAKE_COMMAND} -E copy_directory
+                       ${SRC_DIR} ${BIN_DIR} )
   add_custom_command ( TARGET copy_test_input_files POST_BUILD
                        COMMAND ${CMAKE_COMMAND} -E remove -f
-                       ${BIN_DIR}/*.md ${BIN_DIR}/CMakeLists.*
-                       ${BIN_DIR}/*/*.md ${BIN_DIR}/*/*.png ${BIN_DIR}/*/CMakeLists.* )
+                       ${BIN_DIR}/*.md ${BIN_DIR}/*/*.md
+                       ${BIN_DIR}/CMakeLists.* ${BIN_DIR}/*/CMakeLists.* )
+  if ( "${TEST}" STREQUAL "InversePy" )
+    add_custom_command ( TARGET copy_test_input_files POST_BUILD
+                         COMMAND ${CMAKE_COMMAND} -E remove -f
+                         ${BIN_DIR}/*.py ${BIN_DIR}/*/run_API.py )
+    add_custom_command ( TARGET copy_test_input_files POST_BUILD
+                         COMMAND ${CMAKE_COMMAND} -E remove_directory ${BIN_DIR}/img )
+  endif ( "${TEST}" STREQUAL "InversePy" )
 endforeach ( TEST ${TEST_DIRS} )
diff --git a/InversePy/beam_longitudinal/README.md b/InversePy/beam_longitudinal/README.md
index ff4b299..d750a6b 100644
--- a/InversePy/beam_longitudinal/README.md
+++ b/InversePy/beam_longitudinal/README.md
@@ -40,4 +40,4 @@ The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 The results satisfy the criterion of acceptance.
 
-![Longitudinal beam](longitudinal.PNG "Longitudinal beam response")
+![Longitudinal beam](../img/longitudinal.png "Longitudinal beam response")
diff --git a/InversePy/beam_longitudinal/run_API.py b/InversePy/beam_longitudinal/run_API.py
index 33392b0..c134636 100644
--- a/InversePy/beam_longitudinal/run_API.py
+++ b/InversePy/beam_longitudinal/run_API.py
@@ -11,13 +11,13 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     init_length = 1.0
 
     twin_funcId = [5]  # 5 .. rel distance triad 113-112
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     # holds the results from the base data (graphical/numerical comparison)
     baseRes = []
@@ -31,9 +31,6 @@ def run_simulation(inp):
 
     # ======================= TIME LOOP ==========================
     for n in range(inp["total_increments"]):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
-
         base[0] = baseD[n] - init_length  # data modification based on base output
         baseRes.append([baseD[n]])
 
@@ -43,13 +40,13 @@ def run_simulation(inp):
 
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
-    baseR, twinR = run_simulation(utils.read_input())
+    baseR, twinR = run_simulation(utils.read_input(), print_stepping)
 
     # graphical output with matplotlib
     if plot_curves:
@@ -64,6 +61,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/beam_modal_forces/README.md b/InversePy/beam_modal_forces/README.md
index 0559ff0..06dec80 100644
--- a/InversePy/beam_modal_forces/README.md
+++ b/InversePy/beam_modal_forces/README.md
@@ -37,4 +37,4 @@ The curve plots show the measurements (reference) as line,
 the results from the inverse method are presented as points.
 The results satisfy the criterion of acceptance.
 
-![Transversal beam deflection](modal_forces.png "Transversal beam deflection")
+![Transversal beam deflection](../img/modal_forces.png "Transversal beam deflection")
diff --git a/InversePy/beam_modal_forces/run_API.py b/InversePy/beam_modal_forces/run_API.py
index e9fdbd1..82bb5b1 100644
--- a/InversePy/beam_modal_forces/run_API.py
+++ b/InversePy/beam_modal_forces/run_API.py
@@ -11,11 +11,11 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     twin_funcId = [1, 3, 4]  # displacements in DY_17, DY_113 and DY_114
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     # displacements from file
     baseRes = utils.read_data_from_file("./refData.asc")
@@ -24,8 +24,6 @@ def run_simulation(inp):
 
     # ======================= TIME LOOP ==========================
     for n in range(inp["total_increments"]):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
 
         base = baseRes[n, 0:2]  # prepare input for twin, only internal forces
 
@@ -35,13 +33,13 @@ def run_simulation(inp):
 
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
-    baseR, twinR = run_simulation(utils.read_input())
+    baseR, twinR = run_simulation(utils.read_input(), print_stepping)
 
     # graphical output with matplotlib
     if plot_curves:
@@ -54,6 +52,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/beam_transversal/README.md b/InversePy/beam_transversal/README.md
index b8a9433..171039f 100644
--- a/InversePy/beam_transversal/README.md
+++ b/InversePy/beam_transversal/README.md
@@ -39,6 +39,6 @@ The curve plots show the measurements (reference) as line,
 the results from the inverse method are presented as points.
 The results satisfy the criterion of acceptance.
 
-![Transversal beam moment](transversal_1.PNG "Transversal beam moment")
+![Transversal beam moment](../img/transversal1.png "Transversal beam moment")
 
-![Transversal beam deflection](transversal_2.PNG "Transversal beam deflection")
+![Transversal beam deflection](../img/transversal2.png "Transversal beam deflection")
diff --git a/InversePy/beam_transversal/run_API.py b/InversePy/beam_transversal/run_API.py
index 3d27930..46abe6b 100644
--- a/InversePy/beam_transversal/run_API.py
+++ b/InversePy/beam_transversal/run_API.py
@@ -11,19 +11,17 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     twin_funcId = [8, 9, 4]  # results Mz_17, Qy_113 and Disp_114
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     baseRes = utils.read_data_from_file("./refData.asc", sep=",")
     twinRes = []
 
     # ======================= TIME LOOP ==========================
     for n in range(inp["total_increments"]):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
 
         base = baseRes[n, 0:2]  # prepare input for twin, only internal forces
 
@@ -33,11 +31,11 @@ def run_simulation(inp):
 
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
     # Check if yaml file exists, if not use direct input
     try:
@@ -50,7 +48,7 @@ def main():
         )
         inp['total_increments'] = 1000
 
-    baseR, twinR = run_simulation(inp)
+    baseR, twinR = run_simulation(inp, print_stepping)
 
     # graphical output with matplotlib
     # output internal forces
@@ -73,6 +71,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/corner/README.md b/InversePy/corner/README.md
index 9cc87fe..4f83da1 100644
--- a/InversePy/corner/README.md
+++ b/InversePy/corner/README.md
@@ -54,13 +54,10 @@ the data read-in, by using the utility method sensor_data_matrix().
 The read-in data columns are specified in the yaml file or by using a
 dictionary as parameter in the method `internal_eq`.
 
-The directory contains a .png file showing the match between the base and
-twin simulation.
-
 # Response data
 
 The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 The results satisfy the criterion of acceptance.
 
-![corner](end_triad_position.png "Displacement response at the end triad")
+![corner](../img/end_triad_position.png "Displacement response at the end triad")
diff --git a/InversePy/corner/run_API.py b/InversePy/corner/run_API.py
index 17e0db2..af77a3e 100644
--- a/InversePy/corner/run_API.py
+++ b/InversePy/corner/run_API.py
@@ -15,20 +15,20 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     # id 2 deflection, id 3 displacement end triad
     twin_funcId = [2, 3]
 
     # initialize analysis model
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     twinRes = []
     baseRes = []
 
     # read the column(s) as specified in the input
     # here column 3 (contains the spring force)
-    dataMat = utils.sensor_data_matrix(inp)
+    dataMat = utils.sensor_data_matrix(inp, print_stepping)
 
     # As an alternative - reading from file refData.asc
     # [0] spring deflection,
@@ -45,8 +45,6 @@ def run_simulation(inp):
 
     # ======================= TIME LOOP =========================
     for n in range(len(dataMat)-1):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
 
         # extract data from data matrix dataMat (spring force)
         force_data = dataMat[n+1,:]
@@ -59,12 +57,12 @@ def run_simulation(inp):
 
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseDisp, baseRes, twinRes
 
 
-def main():
-    # Check if yaml file exixts, if not use direct input
+def main(print_stepping=False):
+    # Check if yaml file exists, if not use direct input
     try:
         inp = utils.read_input()
     except:
@@ -75,7 +73,7 @@ def main():
         )
 
     # Run inverse method (static case)
-    bdsp, baseR, twinR = run_simulation(inp)
+    bdsp, baseR, twinR = run_simulation(inp, print_stepping)
 
     # take second column from baseR (contains the displacements of the end triad)
     tdisp = [ [tr[1]] for tr in twinR]
@@ -101,6 +99,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/crane/README.md b/InversePy/crane/README.md
index 0ecf4bb..127da33 100644
--- a/InversePy/crane/README.md
+++ b/InversePy/crane/README.md
@@ -52,9 +52,6 @@ the spring deflection and the spring force in z direction at the
 revolute joint (several data are given for controlling the functionality, column 1 and
 3 must have the same measurements).
 
-In the directory is a .png file showing the match between the base and
-twin simulation.
-
 # Response data
 
 The curve plot shows the measurements (reference) as line,
@@ -67,4 +64,4 @@ The results satisfy the criterion of acceptance.
 
 The figure below shows a section of the analysis.
 
-![crane_detail](results_zoomView.png "Displacement response at the end triad")
+![crane_detail](../img/results_zoomView.png "Displacement response at the end triad")
diff --git a/InversePy/crane/run_API.py b/InversePy/crane/run_API.py
index 84d2c7a..0cd2a17 100644
--- a/InversePy/crane/run_API.py
+++ b/InversePy/crane/run_API.py
@@ -16,12 +16,12 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     twin_funcId = [1]  # revolute joint, spring force
 
     # initialize model
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     # read in data from base run, file contains 4 colums
     # [0] strut_tz, revolute joint, z trans, length
@@ -40,17 +40,17 @@ def run_simulation(inp):
 
     # ======================= TIME LOOP =========================
     for n in range(inp['total_increments']):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment {}".format(n))
 
         # input joint spring force
         base_n[0] = base[n+1][3]
 
-        print("before: ", base_n)
+        if print_stepping:
+            print("before:", base_n)
         # convert revolute Joint spring force to deflection
         jointIDs = twinmodel.convert_rev_joint_force(base_n)
-        print("Convert spring force to displacement for rev joints: ", jointIDs)
-        print("after: ", base_n)
+        if print_stepping:
+            print("Convert spring force to displacement for rev joints:", jointIDs)
+            print("after:", base_n)
 
         baseSpr.append([base_n[0]])
 
@@ -60,11 +60,11 @@ def run_simulation(inp):
 
         twinSpr.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseDisp, baseSpr, twinSpr
 
 
-def main():
+def main(print_stepping=False):
 
     try:
         inp = utils.read_input()
@@ -76,7 +76,7 @@ def main():
         )
         inp['total_increments'] = 10000
 
-    baseDisp, baseR, twinR = run_simulation(inp)
+    baseDisp, baseR, twinR = run_simulation(inp, print_stepping)
 
     # graphical output with matplotlib for forces
     if plot_curves:
@@ -92,4 +92,4 @@ def main():
     print("----\nCrane example successfully finished\n----")
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
diff --git a/InversePy/frame_disp/README.md b/InversePy/frame_disp/README.md
index 2ddb4d5..9624bcd 100644
--- a/InversePy/frame_disp/README.md
+++ b/InversePy/frame_disp/README.md
@@ -54,4 +54,4 @@ available (by setting `plotCurves=True` in the `run_API.py` script).
 The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 
-![Frame displacement](disp.png "Frame tip displacement")
+![Frame displacement](../imp/frame_disp.png "Frame tip displacement")
diff --git a/InversePy/frame_disp/run_API.py b/InversePy/frame_disp/run_API.py
index 65cbb51..a1f9e6e 100644
--- a/InversePy/frame_disp/run_API.py
+++ b/InversePy/frame_disp/run_API.py
@@ -28,11 +28,11 @@ def run_reducer():
     return ierr
 
 
-def run_simulation(inp, useFtnInv=False):
+def run_simulation(inp, print_stepping, useFtnInv=False):
 
     twin_funcId = [4]
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     base1 = [0.0]
     baseD = utils.read_data_from_file("./refData.asc")
@@ -42,8 +42,6 @@ def run_simulation(inp, useFtnInv=False):
 
     # ======================= TIME LOOP =========================
     for n in range(inp["total_increments"]):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
 
         base1[0] = baseD[n] - 0.5
         baseRes.append([baseD[n]])
@@ -57,11 +55,11 @@ def run_simulation(inp, useFtnInv=False):
 
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
     if path.isfile("frame_S.fmx"):
         print("Reduced FE part found - reusing it")
@@ -72,10 +70,10 @@ def main():
 
     if environ["USE_FTNINV"] == "1":
         print("\n>>> Using the inverse method implemented in Fortran <<<\n")
-        baseR, twinR = run_simulation(utils.read_input(), True)
+        baseR, twinR = run_simulation(utils.read_input(), False, True)
     else:
         print("\n>>> Using the inverse method implemented in Python <<<\n")
-        baseR, twinR = run_simulation(utils.read_input())
+        baseR, twinR = run_simulation(utils.read_input(), print_stepping)
 
     # graphical output with matplotlib
     if plot_curves:
@@ -90,6 +88,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/frame_intF/README.md b/InversePy/frame_intF/README.md
index 27642a6..8a024f9 100644
--- a/InversePy/frame_intF/README.md
+++ b/InversePy/frame_intF/README.md
@@ -53,4 +53,4 @@ available (by setting `plotCurves=True` in the `run_API.py` script).
 The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 
-![Frame internal forces](forces.png "Frame sectional forces")
+![Frame internal forces](../img/frame_forces.png "Frame sectional forces")
diff --git a/InversePy/frame_intF/run_API.py b/InversePy/frame_intF/run_API.py
index 74af601..9a7c9b3 100644
--- a/InversePy/frame_intF/run_API.py
+++ b/InversePy/frame_intF/run_API.py
@@ -11,11 +11,11 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     twin_funcId = [6, 5, 7, 4]
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     baseRes = utils.read_data_from_file("./refData.asc")
 
@@ -23,8 +23,6 @@ def run_simulation(inp):
 
     # ======================= TIME LOOP =========================
     for n in range(inp["total_increments"]):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
 
         base1 = baseRes[n + 1]
 
@@ -34,15 +32,15 @@ def run_simulation(inp):
 
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
     inp = utils.read_input()
 
-    baseR, twinR = run_simulation(inp)
+    baseR, twinR = run_simulation(inp, print_stepping)
 
     # graphical output with matplotlib for forces
     if plot_curves:
@@ -57,6 +55,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/frame_modes/README.md b/InversePy/frame_modes/README.md
index 7a6cfb2..7c67a50 100644
--- a/InversePy/frame_modes/README.md
+++ b/InversePy/frame_modes/README.md
@@ -57,4 +57,4 @@ available (by setting `plotCurves=True` in the `run_API.py` script).
 The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 
-![Frame modes](disp.png "Frame tip displacement")
+![Frame modes](../img/frame_modes_disp.png "Frame tip displacement")
diff --git a/InversePy/frame_modes/run_API.py b/InversePy/frame_modes/run_API.py
index 133f92b..ff0b937 100644
--- a/InversePy/frame_modes/run_API.py
+++ b/InversePy/frame_modes/run_API.py
@@ -11,11 +11,11 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     twin_funcId = [4]
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     base1 = [0.0]
     baseD = utils.read_data_from_file("./refData.asc")
@@ -25,8 +25,6 @@ def run_simulation(inp):
 
     # ======================= TIME LOOP =========================
     for n in range(inp["total_increments"]):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
 
         base1[0] = baseD[n] - 0.5
         baseRes.append([baseD[n]])
@@ -37,17 +35,17 @@ def run_simulation(inp):
 
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
     if not have_sci_py:
         print("----\n FRAME - inverse modes, skipped due to absent scipy")
         return None
 
-    baseR, twinR = run_simulation(utils.read_input())
+    baseR, twinR = run_simulation(utils.read_input(), print_stepping)
 
     # graphical output with matplotlib
     if plot_curves:
@@ -62,6 +60,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/frame_relD/README.md b/InversePy/frame_relD/README.md
index 1f1165b..7b79a83 100644
--- a/InversePy/frame_relD/README.md
+++ b/InversePy/frame_relD/README.md
@@ -53,4 +53,4 @@ available (by setting `plotCurves=True` in the `run_API.py` script).
 The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 
-![Frame relative distance](dispForce.png "Frame tip displacement")
+![Frame relative distance](../img/dispForce.png "Frame tip displacement")
diff --git a/InversePy/frame_relD/run_API.py b/InversePy/frame_relD/run_API.py
index 0039652..8a325e5 100644
--- a/InversePy/frame_relD/run_API.py
+++ b/InversePy/frame_relD/run_API.py
@@ -11,14 +11,14 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     L_0 = 0.496158715  # Initial length
     X_0 = 0.3  # initial position
 
     twin_funcId = [8, 9]
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     base1 = [0.0, 0.0]
 
@@ -28,8 +28,6 @@ def run_simulation(inp):
 
     # ======================= TIME LOOP =========================
     for n in range(inp["total_increments"]):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
 
         base1[:] = baseRes[n + 1]
         base1[0] -= L_0  # subtract initial length
@@ -41,13 +39,13 @@ def run_simulation(inp):
 
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
-    baseR, twinR = run_simulation(utils.read_input())
+    baseR, twinR = run_simulation(utils.read_input(), print_stepping)
 
     # graphical output with matplotlib
     if plot_curves:
@@ -63,6 +61,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/frame_spring/README.md b/InversePy/frame_spring/README.md
index 8384941..f4f63cf 100644
--- a/InversePy/frame_spring/README.md
+++ b/InversePy/frame_spring/README.md
@@ -54,4 +54,4 @@ available (by setting `plot_curves=True` in the `run_API.py` script).
 The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 
-![Frame spring deflection](spring_defl_base-twin.png "Frame tip displacement")
+![Frame spring deflection](../img/spring_defl_base-twin.png "Frame tip displacement")
diff --git a/InversePy/frame_spring/run_API.py b/InversePy/frame_spring/run_API.py
index ead56f8..527c0c4 100644
--- a/InversePy/frame_spring/run_API.py
+++ b/InversePy/frame_spring/run_API.py
@@ -11,12 +11,12 @@
 plot_curves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     # 8 rel Distance, 9 displacement in y and 12 spring deflection
     twin_funcId = [8,9,12]
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     # readin the results from the base simulation (or senor values)
     baseRes = utils.read_data_from_file("./refData.asc")
@@ -27,27 +27,27 @@ def run_simulation(inp):
 
     # ======================= TIME LOOP =========================
     for n in range(len(baseRes)-1):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment {}".format(n))
 
         # base contains the second and third column from the input
         # 2nd colum disp y at the end triad
         # 3rd column spring deflection
         base = [ baseRes[n+1][1], baseRes[n+1][2] ]
         base[0] -= 0.3  # subtract initial position
-        print("Updated measurements: ", base)
+        if print_stepping:
+            print("Updated measurements:", base)
 
         # Inverse method
         twin = twinmodel.run_inverse(base,out_def=twin_funcId)
         if twin is None:
             break; # end of simulation
+
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
     # Check if yaml file exists, if not use direct input
     try:
         inp = utils.read_input()
@@ -59,7 +59,7 @@ def main():
             },
         )
 
-    baseR, twinR = run_simulation(inp)
+    baseR, twinR = run_simulation(inp, print_stepping)
 
     # graphical output with matplotlib
     if plot_curves:
@@ -75,6 +75,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/frame_strain/README.md b/InversePy/frame_strain/README.md
index 67c7150..6e8a25d 100644
--- a/InversePy/frame_strain/README.md
+++ b/InversePy/frame_strain/README.md
@@ -53,4 +53,4 @@ available (by setting `plotCurves=True` in the `run_API.py` script).
 The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 
-![Frame strain](disp.png "Frame tip displacement")
+![Frame strain](../img/disp.png "Frame tip displacement")
diff --git a/InversePy/frame_strain/run_API.py b/InversePy/frame_strain/run_API.py
index 7c0bdef..4ac903c 100644
--- a/InversePy/frame_strain/run_API.py
+++ b/InversePy/frame_strain/run_API.py
@@ -11,19 +11,17 @@
 plotCurves = False  # False for no plots
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     twin_funcId = [10, 7, 4]  # gage bottom (first tensor component), force X, disp Z
 
-    twinmodel = FedemRun(".", inp)
+    twinmodel = FedemRun(".", inp, print_stepping)
 
     baseRes = utils.read_data_from_file("./refData.asc")
     twinRes = []
 
     # ======================= TIME LOOP =========================
     for n in range(inp["total_increments"]):
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
 
         base1 = baseRes[n]
 
@@ -31,14 +29,15 @@ def run_simulation(inp):
         if twin is None:
             break  # end of simulation
 
-        print("twin sensors:", twin)
+        if print_stepping:
+            print("twin sensors:", twin)
         twinRes.append(twin)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
     try:
         inp = utils.read_input()
@@ -52,7 +51,7 @@ def main():
         )
         inp['total_increments'] = 100
 
-    baseR, twinR = run_simulation(inp)
+    baseR, twinR = run_simulation(inp, print_stepping)
 
     # graphical output with matplotlib
     if plotCurves:
@@ -69,6 +68,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/frame_strain/disp.png b/InversePy/img/disp.png
similarity index 100%
rename from InversePy/frame_strain/disp.png
rename to InversePy/img/disp.png
diff --git a/InversePy/frame_relD/dispForce.png b/InversePy/img/dispForce.png
similarity index 100%
rename from InversePy/frame_relD/dispForce.png
rename to InversePy/img/dispForce.png
diff --git a/InversePy/frame_spring/dispY_base-twin.png b/InversePy/img/dispY_base-twin.png
similarity index 100%
rename from InversePy/frame_spring/dispY_base-twin.png
rename to InversePy/img/dispY_base-twin.png
diff --git a/InversePy/corner/end_triad_position.png b/InversePy/img/end_triad_position.png
similarity index 100%
rename from InversePy/corner/end_triad_position.png
rename to InversePy/img/end_triad_position.png
diff --git a/InversePy/frame_strain/force.png b/InversePy/img/force.png
similarity index 100%
rename from InversePy/frame_strain/force.png
rename to InversePy/img/force.png
diff --git a/InversePy/frame_disp/disp.png b/InversePy/img/frame_disp.png
similarity index 100%
rename from InversePy/frame_disp/disp.png
rename to InversePy/img/frame_disp.png
diff --git a/InversePy/frame_intF/forces.png b/InversePy/img/frame_forces.png
similarity index 100%
rename from InversePy/frame_intF/forces.png
rename to InversePy/img/frame_forces.png
diff --git a/InversePy/frame_modes/disp.png b/InversePy/img/frame_modes_disp.png
similarity index 100%
rename from InversePy/frame_modes/disp.png
rename to InversePy/img/frame_modes_disp.png
diff --git a/InversePy/beam_longitudinal/longitudinal.PNG b/InversePy/img/longitudinal.png
similarity index 100%
rename from InversePy/beam_longitudinal/longitudinal.PNG
rename to InversePy/img/longitudinal.png
diff --git a/InversePy/beam_modal_forces/modal_forces.png b/InversePy/img/modal_forces.png
similarity index 100%
rename from InversePy/beam_modal_forces/modal_forces.png
rename to InversePy/img/modal_forces.png
diff --git a/InversePy/frame_relD/relDist.png b/InversePy/img/relDist.png
similarity index 100%
rename from InversePy/frame_relD/relDist.png
rename to InversePy/img/relDist.png
diff --git a/InversePy/crane/results.png b/InversePy/img/results.png
similarity index 100%
rename from InversePy/crane/results.png
rename to InversePy/img/results.png
diff --git a/InversePy/crane/results_zoomView.png b/InversePy/img/results_zoomView.png
similarity index 100%
rename from InversePy/crane/results_zoomView.png
rename to InversePy/img/results_zoomView.png
diff --git a/InversePy/shell_strain/shell_strain.PNG b/InversePy/img/shell_strain.png
similarity index 100%
rename from InversePy/shell_strain/shell_strain.PNG
rename to InversePy/img/shell_strain.png
diff --git a/InversePy/frame_spring/spring_defl_base-twin.png b/InversePy/img/spring_defl_base-twin.png
similarity index 100%
rename from InversePy/frame_spring/spring_defl_base-twin.png
rename to InversePy/img/spring_defl_base-twin.png
diff --git a/InversePy/frame_strain/strain.png b/InversePy/img/strain.png
similarity index 100%
rename from InversePy/frame_strain/strain.png
rename to InversePy/img/strain.png
diff --git a/InversePy/beam_transversal/transversal_1.PNG b/InversePy/img/transversal1.png
similarity index 100%
rename from InversePy/beam_transversal/transversal_1.PNG
rename to InversePy/img/transversal1.png
diff --git a/InversePy/beam_transversal/transversal_2.PNG b/InversePy/img/transversal2.png
similarity index 100%
rename from InversePy/beam_transversal/transversal_2.PNG
rename to InversePy/img/transversal2.png
diff --git a/InversePy/shell_strain/README.md b/InversePy/shell_strain/README.md
index 3784c86..f414db0 100644
--- a/InversePy/shell_strain/README.md
+++ b/InversePy/shell_strain/README.md
@@ -47,4 +47,4 @@ The curve plot shows the measurements (reference) as line,
 the results from the inverse method are presented as points.
 The results satisfy the criterion of acceptance.
 
-![Shell strain](shell_strain.PNG "Shell deflection")
+![Shell strain](../img/shell_strain.png "Shell deflection")
diff --git a/InversePy/shell_strain/run_API.py b/InversePy/shell_strain/run_API.py
index e801085..08ee869 100644
--- a/InversePy/shell_strain/run_API.py
+++ b/InversePy/shell_strain/run_API.py
@@ -28,14 +28,14 @@ def run_reducer():
     return ierr
 
 
-def run_simulation(inp):
+def run_simulation(inp, print_stepping):
 
     twin_funcId = [3]
 
     twinRes = [[0.0]]
     baseRes = []
 
-    twinmodel = FedemRun("./", inp)
+    twinmodel = FedemRun("./", inp, print_stepping)
 
     m_strain = utils.read_data_from_file("./refStrainData.asc")
     baseDisp = utils.read_data_from_file("./refBaseData.asc")
@@ -43,13 +43,11 @@ def run_simulation(inp):
     # ======================= TIME LOOP ==========================
     for n in range(inp["total_increments"]):
 
-        print("\n+++++++++++++++++++++++++++++++++++++++++++++")
-        print("Solving time increment", n)
-
         baseRes.append([baseDisp[n]])
 
         strains_on_thing = m_strain[n]
-        print("strains_on_thing:", strains_on_thing)
+        if print_stepping:
+            print("strains_on_thing:", strains_on_thing)
 
         twinDisp = twinmodel.run_inverse(strains_on_thing, out_def=twin_funcId)
         if twinDisp is None:
@@ -57,11 +55,11 @@ def run_simulation(inp):
 
         twinRes.append(twinDisp)
 
-    twinmodel.solver_done()
+    twinmodel.done_inverse()
     return baseRes, twinRes
 
 
-def main():
+def main(print_stepping=False):
 
     if path.isfile("shell_S.fmx"):
         print("Reduced FE part found - reusing it")
@@ -70,7 +68,7 @@ def main():
         if ierr != 0:
             exit(ierr)
 
-    baseRes, twinRes = run_simulation(utils.read_input())
+    baseRes, twinRes = run_simulation(utils.read_input(), print_stepping)
 
     # graphical output with matplotlib
     if plotCurves:
@@ -85,6 +83,6 @@ def main():
 
 
 if __name__ == "__main__":
-    main()
+    main(not utils.parse_input().no_print)
 
 # end of file
diff --git a/InversePy/test_utils.py b/InversePy/test_utils.py
index f6660cc..0da98f3 100644
--- a/InversePy/test_utils.py
+++ b/InversePy/test_utils.py
@@ -170,6 +170,12 @@ def parse_input(src_dir=None):
         default=0,
         help="Specify runtype: 0 parallel run, 1 serial run, 2 only base, 3 only twin",
     )
+    parser.add_argument(
+        "-q",
+        "--no-print",
+        action="store_true",
+        help="Suppress time step print.",
+    )
     parser.add_argument("-H", action="store_true", help="Print help for current module")
 
     return parser.parse_args()
@@ -183,7 +189,7 @@ def read_input(src_dir=None):
     return safe_load(open(args.input, "r"))[args.method][args.case]
 
 
-def sensor_data_matrix(inp):  # NOSONAR
+def sensor_data_matrix(inp, do_print=True):  # NOSONAR
     """
     Parsing data source
     Supports only the reading from file!
@@ -200,7 +206,8 @@ def sensor_data_matrix(inp):  # NOSONAR
     ]
     eq_list_def = [key for key in inp["internal_equations"] if key in eq_list]
 
-    print("inp:",inp)
+    if do_print:
+        print("inp:",inp)
 
     data_mat = None
     for item in eq_list_def:
@@ -233,7 +240,7 @@ def sensor_data_matrix(inp):  # NOSONAR
                             raise ValueError("Only reading from file is supported")
                         break
 
-    if data_mat is not None:
+    if data_mat is not None and do_print:
         print("\ntype: ", type(data_mat), " --- size: ", np.size(data_mat, axis = 0))
 
     return data_mat