diff --git a/DG_Approximation.py b/DG_Approximation.py
index cab56cfed73ef70578df114f97a9f152309b291f..2314e6b614666c4e8f37096f13f48ee001dd5dd5 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -3,13 +3,13 @@
 @author: Laura C. Kühle
 
 Urgent:
+TODO: Add way of saving data (np.savez('data/' + name,
+    coefficients=projection, troubled_cells=troubled_cells)) -> Done
+TODO: Move plotting into separate function
 TODO: Adapt TCD from Soraya
     (Dropbox->...->TEST_troubled-cell-detector->Troubled_Cell_Detector)
-TODO: Add way of saving data (np.savez('data/' + name,
-    coefficients=projection, troubled_cells=troubled_cells) )
 TODO: Add verbose output
 TODO: Improve file naming (e.g. use '.' instead of '__')
-TODO: Move plotting into separate function
 
 Critical, but not urgent:
 TODO: Use cfl_number for updating, not just time
@@ -36,6 +36,7 @@ TODO: Add type annotations to function heads
 
 """
 import os
+import json
 import numpy as np
 from sympy import Symbol
 import math
@@ -53,6 +54,18 @@ matplotlib.use('Agg')
 x = Symbol('x')
 
 
+def encode_ndarray(obj):
+    if isinstance(obj, np.ndarray):
+        return obj.tolist()
+    return obj
+
+
+def decode_ndarray(obj):
+    if isinstance(obj, list):
+        return np.asarray(obj)
+    return obj
+
+
 class DGScheme:
     """Class for Discontinuous Galerkin Method.
 
@@ -190,7 +203,7 @@ class DGScheme:
             self._polynomial_degree, self._num_grid_cells, self._detector,
             self._limiter)
 
-    def approximate(self):
+    def approximate(self, data_name):
         """Approximates projection.
 
         Initializes projection and evolves it in time. Each time step consists
@@ -200,6 +213,11 @@ class DGScheme:
         At final time, result and error plots are
         generated and, if verbose flag is set, also displayed.
 
+        Attributes
+        ----------
+        data_name : str
+            Name of data.
+
         """
         projection = self._do_initial_projection(self._init_cond)
 
@@ -227,10 +245,30 @@ class DGScheme:
 
             current_time += time_step
 
+        # Save approximation results in dictionary
+        approx_stats = {'projection': projection, 'time_history': time_history,
+                        'troubled_cell_history': troubled_cell_history}
+
+        # Encode all ndarrays to fit JSON format
+        approx_stats = {key: encode_ndarray(approx_stats[key])
+                        for key in approx_stats.keys()}
+
+        # Save approximation results in JSON format
+        with open(self._plot_dir+'/' + data_name + '.json', 'w') \
+                as json_file:
+            json_file.write(json.dumps(approx_stats))
+
+        # Read approximation results
+        with open(self._plot_dir+'/' + data_name + '.json') as json_file:
+            approx_stats = json.load(json_file)
+
+        # Decode all ndarrays by converting lists
+        approx_stats = {key: decode_ndarray(approx_stats[key])
+                        for key in approx_stats.keys()}
+
         # Plot exact/approximate results, errors, shock tubes,
         # and any detector-dependant plots
-        self._detector.plot_results(projection, troubled_cell_history,
-                                    time_history)
+        self._detector.plot_results(**approx_stats)
 
     def save_plots(self, plot_name):
         """Saves plotted results.
diff --git a/workflows/approximation.smk b/workflows/approximation.smk
index 293d0f48752d54ae0fa280f6702efe3202acd83c..fa509bbc75d103943d41d1a352b88e722af0d2fa 100644
--- a/workflows/approximation.smk
+++ b/workflows/approximation.smk
@@ -48,7 +48,7 @@ rule approximate_solution:
             print(params.dg_params)
             dg_scheme = DGScheme(plot_dir=params.plot_dir, **params.dg_params)
 
-            dg_scheme.approximate()
+            dg_scheme.approximate(wildcards.scheme)
             dg_scheme.save_plots(wildcards.scheme)
 
             toc = time.perf_counter()