diff --git a/dataprocessing/readtools.py b/dataprocessing/readtools.py
index 91df32902a40416683192dade1e42a12fa3d03dc..3b302feb025a95aa0d6ae51e3a1929d3ad68bcad 100644
--- a/dataprocessing/readtools.py
+++ b/dataprocessing/readtools.py
@@ -1,14 +1,25 @@
 import numpy as np
 import pandas as pd
 from .timeseries import *
+import re
 
 
-def read_timeseries_Modelica(filename, timeseries_names=None):
+def read_timeseries_Modelica(filename, timeseries_names=None, is_regex=False):
     from modelicares import SimRes
     sim = SimRes(filename)
-    if timeseries_names is None:
-        # No trajectory names specified, thus read in all
-        print('TBD')
+    if timeseries_names is None and is_regex is False:
+        # No trajectory names or regex specified, thus read in all
+        timeseries = []
+        for name in sim.names():
+            timeseries.append(TimeSeries(name, sim(name).times(), sim(name).values()))
+    elif is_regex is True:
+        # Read in variables which match with regex
+        timeseries = []
+        p = re.compile(timeseries_names)
+        timeseries_names = [name for name in sim.names() if p.search(name)]
+        timeseries_names.sort()
+        for name in timeseries_names:
+            timeseries.append(TimeSeries(name, sim(name).times(), sim(name).values()))
     else:
         # Read in specified time series
         if not isinstance(timeseries_names, list):
diff --git a/dataprocessing/timeseries.py b/dataprocessing/timeseries.py
index 7115f4ccfaa9bd7964f800cb35c4b3f7414e40a3..83b979ac8c1001a5acef7ddee943ba98a443c7f6 100644
--- a/dataprocessing/timeseries.py
+++ b/dataprocessing/timeseries.py
@@ -14,9 +14,14 @@ class TimeSeries:
     @staticmethod
     def diff(name, ts1, ts2):
         """Returns difference between values of two Timeseries objects.
-        Assumes the same time steps for both timeseries.
         """
-        ts_diff = TimeSeries(name, ts1.time, (ts1.values - ts2.values))
+        if ts1.time==ts2.time:
+            ts_diff = TimeSeries(name, ts1.time, (ts1.values - ts2.values))
+        else:  # different timestamps, common time vector and interpolation required before substraction
+            time = sorted(set(list(ts1.time) + list(ts2.time)))
+            interp_vals_ts1 = np.interp(time, ts1.time, ts1.values)
+            interp_vals_ts2 = np.interp(time, ts2.time, ts2.values)
+            ts_diff = TimeSeries(name, time, (interp_vals_ts2 - interp_vals_ts1))
         return ts_diff
 
 
diff --git a/examples/ExamplesReadModelica.py b/examples/ExamplesReadModelica.py
new file mode 100644
index 0000000000000000000000000000000000000000..eab9a32d6f3292e9a4b1c96362c87a97c5cd1636
--- /dev/null
+++ b/examples/ExamplesReadModelica.py
@@ -0,0 +1,45 @@
+from dataprocessing.readtools import *
+from dataprocessing.plottools import *
+import matplotlib.pyplot as plt
+
+
+# Example 1: read in single variable included in the Modelica results file
+voltage_node126 = read_timeseries_Modelica(
+    r"\\tsclient\N\Research\German Public\ACS0049_SINERGIEN_bsc\Data\WorkData\SimulationResults\IEEE European\Single_scenario_fixed_PV\IEEEEuropean_60.mat",
+    timeseries_names="N126.Vrel")
+plt.figure(1, figsize=(12,8))
+set_timeseries_labels(voltage_node126, "voltage N126")
+plt.plot(voltage_node126.time/3600, voltage_node126.values, label=voltage_node126.label)
+plt.legend()
+plt.show(block=True)
+
+# Example 2: read in multiple variables defined in a list
+voltage_two_nodes = read_timeseries_Modelica(
+    r"\\tsclient\N\Research\German Public\ACS0049_SINERGIEN_bsc\Data\WorkData\SimulationResults\IEEE European\Single_scenario_fixed_PV\IEEEEuropean_60.mat",
+    timeseries_names=["N127.Vrel", "N128.Vrel"])
+plt.figure(2, figsize=(12,8))
+plt.plot(voltage_two_nodes[0].time/3600, voltage_two_nodes[0].values, label=voltage_two_nodes[0].label)
+plt.plot(voltage_two_nodes[1].time/3600, voltage_two_nodes[1].values, label=voltage_two_nodes[1].label)
+plt.legend()
+plt.show(block=True)
+
+# Example 3: read in all voltages using regular expressions
+voltages_all_nodes = read_timeseries_Modelica(
+    r"\\tsclient\N\Research\German Public\ACS0049_SINERGIEN_bsc\Data\WorkData\SimulationResults\IEEE European\Single_scenario_fixed_PV\IEEEEuropean_60.mat",
+    timeseries_names='^[^.]*.Vrel$', is_regex=True)
+plt.figure(3, figsize=(12, 8))
+for i in range(len(voltages_all_nodes)):
+    plt.plot(voltages_all_nodes[i].time / 3600, voltages_all_nodes[i].values, label=voltages_all_nodes[i].label)
+plt.legend()
+plt.show(block=True)
+
+# Example 4: read in all variables
+variables_all = read_timeseries_Modelica(
+    r"\\tsclient\N\Research\German Public\ACS0049_SINERGIEN_bsc\Data\WorkData\SimulationResults\IEEE European\Single_scenario_fixed_PV\IEEEEuropean_60.mat")
+dict_variables_all = {}
+for ts in variables_all:
+    dict_variables_all[ts.name] = ts
+plt.figure(4, figsize=(12, 8))
+plt.plot(dict_variables_all["L12.Irel"].time/3600, dict_variables_all["L12.Irel"].values, label=dict_variables_all["L12.Irel"].label)
+plt.legend()
+plt.show(block=True)
\ No newline at end of file
diff --git a/examples/compare_modelica_dpsim.py b/examples/compare_modelica_dpsim.py
index f87bfca1f1ea8b6f336f168bc24b0385a7428de6..28600a13b4a46e5ff8d66694e4a1d57f3edca587 100644
--- a/examples/compare_modelica_dpsim.py
+++ b/examples/compare_modelica_dpsim.py
@@ -3,26 +3,35 @@ from dataprocessing.plottools import *
 import matplotlib.pyplot as plt
 from plottingtools.config import *
 
-current_EMT = read_timeseries_Modelica(    r"C:\Workspace\ReferenceExamples\Modelica\Synchronous Generator\UnitTest_Eremia_3rdOrderModel_Euler_1ms.mat", ["synchronousGenerator_Park.i[1]"])
+current_emt_mod = read_timeseries_Modelica(r"\\tsclient\N\Research\German Public\ACS0049_SINERGIEN_bsc\Data\WorkData\SimulationResults\SynchronousGenerator\EMT\UnitTest_Kundur_FullModel_Euler_1us.mat", ["synchronousGenerator_Park.i[1]"]) # Note: both results include only one damper winding in q axis
+current_emt_dpsim = read_timeseries_dpsim_real(r"\\tsclient\N\Research\German Public\ACS0049_SINERGIEN_bsc\Data\WorkData\SimulationResults\SynchronousGenerator\EMT\DPsim\UnitTest_FullModel_Trap_1us\data_j.csv")[0]
+current_emt_dpsim.values = -current_emt_dpsim.values
 
+# Comparison EMT
 figure_id = 1
 plt.figure(figure_id, figsize=(12,8))
-set_timeseries_labels(current_EMT, ["EMT"])
-plot_timeseries(figure_id, current_EMT, plt_color=blue)
+set_timeseries_labels(current_emt_mod, ["EMT Modelica"])
+plot_timeseries(figure_id, current_emt_mod)
+set_timeseries_labels(current_emt_dpsim, "EMT DPsim") # TODO: modelica timeseries needs list element, dpsim timeseries needs string
+plot_timeseries(figure_id, current_emt_dpsim, plt_linestyle=':')
 plt.xlabel('Zeit [s]')
 plt.ylabel('Strom [A]')
-plt.show(block=True)
+plt.show(block=False)
 
-multi_current_EMT = read_timeseries_Modelica(r"C:\Workspace\ReferenceExamples\Modelica\Synchronous Generator\UnitTest_Eremia_3rdOrderModel_Euler_1ms.mat",[["synchronousGenerator_Park.i[1]"],["synchronousGenerator_Park.i[2]"]])
+# Comparison DP
+current_dp_mod = read_timeseries_Modelica(r"\\tsclient\N\Research\German Public\ACS0049_SINERGIEN_bsc\Data\WorkData\SimulationResults\SynchronousGenerator\DP\UnitTest_Kundur_FullModel_Euler_1us.mat", ["synchronousGenerator_Park.I[1]"]) # Note: both results include only one damper winding in q axis
+current_dp_dpsim = read_timeseries_dpsim_cmpl(r"\\tsclient\N\Research\German Public\ACS0049_SINERGIEN_bsc\Data\WorkData\SimulationResults\SynchronousGenerator\DP\DPsim\UnitTest_FullModel_Trap_1us\data_j.csv")[0]
+current_dp_dpsim.values = -current_dp_dpsim.values
+current_dpabs_dpsim = current_dp_dpsim.abs(current_dp_dpsim.name + ' abs')
 
 figure_id = 2
 plt.figure(figure_id, figsize=(12,8))
-set_timeseries_labels(multi_current_EMT, ["Phase a","Phase b","Phase c"])
-plot_timeseries(figure_id, multi_current_EMT)
+set_timeseries_labels(current_dp_mod, ["DP Modelica"])
+plot_timeseries(figure_id, current_dp_mod)
+set_timeseries_labels(current_dpabs_dpsim, "DP DPsim") # TODO: modelica timeseries needs list element, dpsim timeseries needs string
+plot_timeseries(figure_id, current_dpabs_dpsim, plt_linestyle=':')
 plt.xlabel('Zeit [s]')
 plt.ylabel('Strom [A]')
 plt.show(block=True)
 
 
-
-