diff --git a/.gitignore b/.gitignore
index 5798e24..6bd035b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -23,4 +23,11 @@ docs/site/
 # environment.
 Manifest.toml
 
-.vscode
\ No newline at end of file
+.vscode
+
+#Julia system image
+JuliaSysimage.dll
+
+#python stuff
+__pycache__
+*.pyc
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 206dd51..ce81cc2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.1.2"
 
 [deps]
 ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -15,3 +16,4 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
diff --git a/Python Burgers example/__pycache__/benchmarkedBurgers.cpython-310.pyc b/Python Burgers example/__pycache__/benchmarkedBurgers.cpython-310.pyc
new file mode 100644
index 0000000..72abefb
Binary files /dev/null and b/Python Burgers example/__pycache__/benchmarkedBurgers.cpython-310.pyc differ
diff --git a/Python Burgers example/benchmarkedBurgers.py b/Python Burgers example/benchmarkedBurgers.py
new file mode 100644
index 0000000..e39f88d
--- /dev/null
+++ b/Python Burgers example/benchmarkedBurgers.py	
@@ -0,0 +1,352 @@
+# 2nd-order accurate finite-volume implementation of the inviscid Burger's
+# equation with piecewise linear slope reconstruction
+#
+# We are solving u_t + u u_x = 0 with outflow boundary conditions
+#
+# M. Zingale (2013-03-26)
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import sys
+import time
+import cProfile
+from memory_profiler import profile
+from pympler import asizeof
+from pympler import muppy
+from pympler import summary
+from pympler import classtracker
+import tracemalloc
+import linecache
+import os
+
+mpl.rcParams['mathtext.fontset'] = 'cm'
+mpl.rcParams['mathtext.rm'] = 'serif'
+
+mpl.rcParams['font.size'] = 12
+mpl.rcParams['legend.fontsize'] = 'large'
+mpl.rcParams['figure.titlesize'] = 'medium'
+
+class Grid1d(object):
+    def __init__(self, nx, ng, xmin=0.0, xmax=1.0, bc="outflow"):
+
+        self.nx = nx
+        self.ng = ng
+
+        self.xmin = xmin
+        self.xmax = xmax
+
+        self.bc=bc
+
+        # python is zero-based.  Make easy intergers to know where the
+        # real data lives
+        self.ilo = ng
+        self.ihi = ng+nx-1
+
+        # physical coords -- cell-centered, left and right edges
+        self.dx = (xmax - xmin)/(nx)
+        self.x = xmin + (np.arange(nx+2*ng)-ng+0.5)*self.dx
+
+        # storage for the solution
+        self.u = np.zeros((nx+2*ng), dtype=np.float64)
+
+    def scratch_array(self):
+        """ return a scratch array dimensioned for our grid """
+        return np.zeros((self.nx+2*self.ng), dtype=np.float64)
+
+    def fill_BCs(self):
+        """ fill all ghostcells as periodic """
+
+        if self.bc == "periodic":
+
+            # left boundary
+            self.u[0:self.ilo] = self.u[self.ihi-self.ng+1:self.ihi+1]
+
+            # right boundary
+            self.u[self.ihi+1:] = self.u[self.ilo:self.ilo+self.ng]
+
+        elif self.bc == "outflow":
+
+            # left boundary
+            self.u[0:self.ilo] = self.u[self.ilo]
+
+            # right boundary
+            self.u[self.ihi+1:] = self.u[self.ihi]
+
+        else:
+            sys.exit("invalid BC")
+
+
+    def norm(self, e):
+        """ return the norm of quantity e which lives on the grid """
+        if len(e) != 2*self.ng + self.nx:
+            return None
+
+        return np.sqrt(self.dx*np.sum(e[self.ilo:self.ihi+1]**2))
+
+
+class Simulation(object):
+
+    def __init__(self, grid):
+        self.grid = grid
+        self.t = 0.0
+
+
+    def init_cond(self, type="tophat"):
+
+        if type == "tophat":
+            self.grid.u[np.logical_and(self.grid.x >= 0.333,
+                                          self.grid.x <= 0.666)] = 1.0
+
+        elif type == "sine":
+            self.grid.u[:] = 1.0
+
+            index = np.logical_and(self.grid.x >= 0.333,
+                                      self.grid.x <= 0.666)
+            self.grid.u[index] += \
+                0.5*np.sin(2.0*np.pi*(self.grid.x[index]-0.333)/0.333)
+
+        elif type == "rarefaction":
+            self.grid.u[:] = 1.0
+            self.grid.u[self.grid.x > 0.5] = 2.0
+
+
+    def timestep(self, C):
+        return C*self.grid.dx/max(abs(self.grid.u[self.grid.ilo:
+                                                  self.grid.ihi+1]))
+
+    def states(self, dt):
+        """ compute the left and right interface states """
+
+        g = self.grid
+        # compute the piecewise linear slopes -- 2nd order MC limiter
+        # we pick a range of cells that includes 1 ghost cell on either
+        # side
+        ib = g.ilo-1
+        ie = g.ihi+1
+
+        u = g.u
+
+        # this is the MC limiter from van Leer (1977), as given in
+        # LeVeque (2002).  Note that this is slightly different than
+        # the expression from Colella (1990)
+
+        dc = g.scratch_array()
+        dl = g.scratch_array()
+        dr = g.scratch_array()
+
+        dc[ib:ie+1] = 0.5*(u[ib+1:ie+2] - u[ib-1:ie  ])
+        dl[ib:ie+1] = u[ib+1:ie+2] - u[ib  :ie+1]
+        dr[ib:ie+1] = u[ib  :ie+1] - u[ib-1:ie  ]
+
+        # these where's do a minmod()
+        d1 = 2.0*np.where(np.fabs(dl) < np.fabs(dr), dl, dr)
+        d2 = np.where(np.fabs(dc) < np.fabs(d1), dc, d1)
+        ldeltau = np.where(dl*dr > 0.0, d2, 0.0)
+
+        # now the interface states.  Note that there are 1 more interfaces
+        # than zones
+        ul = g.scratch_array()
+        ur = g.scratch_array()
+
+        # are these indices right?
+        #
+        #  --+-----------------+------------------+
+        #     ^       i       ^ ^        i+1
+        #     ur(i)     ul(i+1) ur(i+1)
+        #
+        ur[ib:ie+2] = u[ib:ie+2] - \
+                      0.5*(1.0 + u[ib:ie+2]*dt/self.grid.dx)*ldeltau[ib:ie+2]
+
+        ul[ib+1:ie+2] = u[ib:ie+1] + \
+                        0.5*(1.0 - u[ib:ie+1]*dt/self.grid.dx)*ldeltau[ib:ie+1]
+
+        return ul, ur
+
+
+    def riemann(self, ul, ur):
+        """
+        Riemann problem for Burgers' equation.
+        """
+
+        S = 0.5*(ul + ur)
+        ushock = np.where(S > 0.0, ul, ur)
+        ushock = np.where(S == 0.0, 0.0, ushock)
+
+        # rarefaction solution
+        urare = np.where(ur <= 0.0, ur, 0.0)
+        urare = np.where(ul >= 0.0, ul, urare)
+
+        us = np.where(ul > ur, ushock, urare)
+
+        return 0.5*us*us
+
+    def update(self, dt, flux):
+        """ conservative update """
+
+        g = self.grid
+
+        unew = g.scratch_array()
+
+        unew[g.ilo:g.ihi+1] = g.u[g.ilo:g.ihi+1] + \
+            dt/g.dx * (flux[g.ilo:g.ihi+1] - flux[g.ilo+1:g.ihi+2])
+
+        return unew
+
+    def evolve(self, C, tmax):
+
+        self.t = 0.0
+
+        g = self.grid
+
+        # main evolution loop
+        while (self.t < tmax):
+
+            # fill the boundary conditions
+            g.fill_BCs()
+
+            # get the timestep
+            dt = self.timestep(C)
+
+            if (self.t + dt > tmax):
+                dt = tmax - self.t
+
+            # get the interface states
+            ul, ur = self.states(dt)
+
+            # solve the Riemann problem at all interfaces
+            flux = self.riemann(ul, ur)
+
+            # do the conservative update
+            unew = self.update(dt, flux)
+
+            self.grid.u[:] = unew[:]
+
+            self.t += dt
+
+def display_top(snapshot, key_type='lineno', limit=10):
+    snapshot = snapshot.filter_traces((
+        tracemalloc.Filter(False, "<frozen importlib._bootstrap>"),
+        tracemalloc.Filter(False, "<unknown>"),
+    ))
+    top_stats = snapshot.statistics(key_type)
+
+    other = top_stats[limit:]
+    if other:
+        size = sum(stat.size for stat in other)
+        print("%s other: %.1f KiB" % (len(other), size / 1024))
+    total = sum(stat.size for stat in top_stats)
+    print("Total allocated size: %.1f KiB" % (total / 1024))
+            
+def main():
+    #-----------------------------------------------------------------------------
+    # sine
+
+    tracemalloc.start()
+    
+    xmin = 0.0
+    xmax = 1.0
+    nx = 256
+    ng = 2
+    g = Grid1d(nx, ng, bc="periodic")
+
+    # maximum evolution time based on period for unit velocity
+    tmax = (xmax - xmin)/1.0
+
+    C = 0.8
+
+    plt.clf()
+    
+    s = Simulation(g)
+
+    start = time.perf_counter()
+
+    for i in range(0, 10):
+        tend = (i+1)*0.02*tmax
+        s.init_cond("sine")
+
+        uinit = s.grid.u.copy()
+
+        s.evolve(C, tend)
+
+        c = 1.0 - (0.1 + i*0.1)
+        g = s.grid
+        #plt.plot(g.x[g.ilo:g.ihi+1], g.u[g.ilo:g.ihi+1], color=str(c))
+
+    timeTaken = time.perf_counter() - start
+    
+    g = s.grid
+    
+    plt.plot(g.x[g.ilo:g.ihi+1], uinit[g.ilo:g.ihi+1], ls=":", color="0.9", zorder=-1)
+
+    xs = g.x[g.ilo:g.ihi+1]
+    #us = uinit[g.ilo:g.ihi+1]
+    us = g.u[g.ilo:g.ihi+1]
+    
+    plt.xlabel("$x$")
+    plt.ylabel("$u$")
+    plt.savefig("fv-burger-sine.pdf")
+
+    f = open('python-xs.txt','w')
+    for i in range(len(xs)):
+        f.write(str(xs[i]) + '\n')
+    f.close()
+
+    f = open('python-us.txt','w')
+    for i in range(len(us)):
+        f.write(str(us[i]) + '\n')
+    f.close()
+
+    snapshot = tracemalloc.take_snapshot()
+
+    display_top(snapshot)
+
+    print(asizeof.asizeof(g))
+    print(asizeof.asizeof(s))
+
+    tracemalloc.stop()
+
+    return timeTaken
+
+
+    #-----------------------------------------------------------------------------
+    # rarefaction
+    """
+    xmin = 0.0
+    xmax = 1.0
+    nx = 256
+    ng = 2
+    g = Grid1d(nx, ng, bc="outflow")
+
+    # maximum evolution time based on period for unit velocity
+    tmax = (xmax - xmin)/1.0
+
+    C = 0.8
+
+    plt.clf()
+
+    s = Simulation(g)
+
+    for i in range(0, 10):
+        tend = (i+1)*0.02*tmax
+
+        s.init_cond("rarefaction")
+
+        uinit = s.grid.u.copy()
+
+        s.evolve(C, tend)
+
+        c = 1.0 - (0.1 + i*0.1)
+        plt.plot(g.x[g.ilo:g.ihi+1], g.u[g.ilo:g.ihi+1], color=str(c))
+
+
+    plt.plot(g.x[g.ilo:g.ihi+1], uinit[g.ilo:g.ihi+1], ls=":", color="0.9", zorder=-1)
+
+    plt.xlabel("$x$")
+    plt.ylabel("$u$")
+
+    plt.savefig("fv-burger-rarefaction.pdf")"""
+
+if __name__ == "__main__":
+    main()
+    
diff --git a/Python Burgers example/fv-burger-sine.pdf b/Python Burgers example/fv-burger-sine.pdf
new file mode 100644
index 0000000..65d48c1
Binary files /dev/null and b/Python Burgers example/fv-burger-sine.pdf differ
diff --git a/Python Burgers example/mprof0 b/Python Burgers example/mprof0
new file mode 100644
index 0000000..87007ae
--- /dev/null
+++ b/Python Burgers example/mprof0	
@@ -0,0 +1,10 @@
+Filename: C:\Users\alice\AppData\Local\Programs\Python\Python310\lib\site-packages\memory_profiler.py
+
+Line #    Mem usage    Increment  Occurrences   Line Contents
+=============================================================
+  1185    132.2 MiB    132.2 MiB           1               @wraps(wrapped=func)
+  1186                                                     def wrapper(*args, **kwargs):
+  1187    132.2 MiB      0.0 MiB           1                   prof = get_prof()
+  1188    132.2 MiB     -0.0 MiB           1                   val = prof(func)(*args, **kwargs)
+  1189    132.2 MiB      0.0 MiB           1                   show_results_bound(prof)
+  1190    132.2 MiB      0.0 MiB           1                   return val
\ No newline at end of file
diff --git a/Python Burgers example/mprofile_20230427223539.dat b/Python Burgers example/mprofile_20230427223539.dat
new file mode 100644
index 0000000..8771864
--- /dev/null
+++ b/Python Burgers example/mprofile_20230427223539.dat	
@@ -0,0 +1,4 @@
+CMDLINE C:\Users\alice\AppData\Local\Programs\Python\Python310\python.exe run-burgers.py
+MEM 1.593750 1682631339.4155
+MEM 12.250000 1682631339.5263
+MEM 18.988281 1682631339.6373
diff --git a/Python Burgers example/mprofile_20230427223549.dat b/Python Burgers example/mprofile_20230427223549.dat
new file mode 100644
index 0000000..5d5cd4c
--- /dev/null
+++ b/Python Burgers example/mprofile_20230427223549.dat	
@@ -0,0 +1,17 @@
+CMDLINE C:\Users\alice\AppData\Local\Programs\Python\Python310\python.exe run-burgers.py
+MEM 1.593750 1682631349.9088
+MEM 12.308594 1682631350.0204
+MEM 18.761719 1682631350.1310
+MEM 26.718750 1682631350.2426
+MEM 33.562500 1682631350.3522
+MEM 41.679688 1682631350.4640
+MEM 47.640625 1682631350.5757
+MEM 51.242188 1682631350.6864
+MEM 54.992188 1682631350.7981
+MEM 60.523438 1682631350.9096
+MEM 61.910156 1682631351.0213
+MEM 68.863281 1682631351.1321
+MEM 73.292969 1682631351.2437
+MEM 77.042969 1682631351.3551
+MEM 77.992188 1682631351.4668
+MEM 76.988281 1682631351.5779
diff --git a/Python Burgers example/mprofile_20230427223610.dat b/Python Burgers example/mprofile_20230427223610.dat
new file mode 100644
index 0000000..34422b7
--- /dev/null
+++ b/Python Burgers example/mprofile_20230427223610.dat	
@@ -0,0 +1,17 @@
+CMDLINE C:\Users\alice\AppData\Local\Programs\Python\Python310\python.exe benchmarkedBurgers.py
+MEM 1.593750 1682631370.7644
+MEM 12.574219 1682631370.8785
+MEM 20.113281 1682631370.9894
+MEM 26.824219 1682631371.1002
+MEM 33.320312 1682631371.2108
+MEM 41.953125 1682631371.3211
+MEM 47.867188 1682631371.4328
+MEM 51.718750 1682631371.5430
+MEM 55.390625 1682631371.6541
+MEM 60.632812 1682631371.7637
+MEM 62.015625 1682631371.8754
+MEM 68.910156 1682631371.9862
+MEM 73.503906 1682631372.0975
+MEM 76.750000 1682631372.2087
+MEM 77.832031 1682631372.3204
+MEM 77.070312 1682631372.4294
diff --git a/Python Burgers example/originalBurgers.py b/Python Burgers example/originalBurgers.py
new file mode 100644
index 0000000..9494fa9
--- /dev/null
+++ b/Python Burgers example/originalBurgers.py	
@@ -0,0 +1,304 @@
+# 2nd-order accurate finite-volume implementation of the inviscid Burger's
+# equation with piecewise linear slope reconstruction
+#
+# We are solving u_t + u u_x = 0 with outflow boundary conditions
+#
+# M. Zingale (2013-03-26)
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import sys
+
+mpl.rcParams['mathtext.fontset'] = 'cm'
+mpl.rcParams['mathtext.rm'] = 'serif'
+
+mpl.rcParams['font.size'] = 12
+mpl.rcParams['legend.fontsize'] = 'large'
+mpl.rcParams['figure.titlesize'] = 'medium'
+
+class Grid1d(object):
+    def __init__(self, nx, ng, xmin=0.0, xmax=1.0, bc="outflow"):
+
+        self.nx = nx
+        self.ng = ng
+
+        self.xmin = xmin
+        self.xmax = xmax
+
+        self.bc=bc
+
+        # python is zero-based.  Make easy intergers to know where the
+        # real data lives
+        self.ilo = ng
+        self.ihi = ng+nx-1
+
+        # physical coords -- cell-centered, left and right edges
+        self.dx = (xmax - xmin)/(nx)
+        self.x = xmin + (np.arange(nx+2*ng)-ng+0.5)*self.dx
+
+        # storage for the solution
+        self.u = np.zeros((nx+2*ng), dtype=np.float64)
+
+
+    def scratch_array(self):
+        """ return a scratch array dimensioned for our grid """
+        return np.zeros((self.nx+2*self.ng), dtype=np.float64)
+
+
+    def fill_BCs(self):
+        """ fill all ghostcells as periodic """
+
+        if self.bc == "periodic":
+
+            # left boundary
+            self.u[0:self.ilo] = self.u[self.ihi-self.ng+1:self.ihi+1]
+
+            # right boundary
+            self.u[self.ihi+1:] = self.u[self.ilo:self.ilo+self.ng]
+
+        elif self.bc == "outflow":
+
+            # left boundary
+            self.u[0:self.ilo] = self.u[self.ilo]
+
+            # right boundary
+            self.u[self.ihi+1:] = self.u[self.ihi]
+
+        else:
+            sys.exit("invalid BC")
+
+
+    def norm(self, e):
+        """ return the norm of quantity e which lives on the grid """
+        if len(e) != 2*self.ng + self.nx:
+            return None
+
+        return np.sqrt(self.dx*np.sum(e[self.ilo:self.ihi+1]**2))
+
+
+class Simulation(object):
+
+    def __init__(self, grid):
+        self.grid = grid
+        self.t = 0.0
+
+
+    def init_cond(self, type="tophat"):
+
+        if type == "tophat":
+            self.grid.u[np.logical_and(self.grid.x >= 0.333,
+                                          self.grid.x <= 0.666)] = 1.0
+
+        elif type == "sine":
+            self.grid.u[:] = 1.0
+
+            index = np.logical_and(self.grid.x >= 0.333,
+                                      self.grid.x <= 0.666)
+            self.grid.u[index] += \
+                0.5*np.sin(2.0*np.pi*(self.grid.x[index]-0.333)/0.333)
+
+        elif type == "rarefaction":
+            self.grid.u[:] = 1.0
+            self.grid.u[self.grid.x > 0.5] = 2.0
+
+
+
+    def timestep(self, C):
+        return C*self.grid.dx/max(abs(self.grid.u[self.grid.ilo:
+                                                  self.grid.ihi+1]))
+
+
+    def states(self, dt):
+        """ compute the left and right interface states """
+
+        g = self.grid
+        # compute the piecewise linear slopes -- 2nd order MC limiter
+        # we pick a range of cells that includes 1 ghost cell on either
+        # side
+        ib = g.ilo-1
+        ie = g.ihi+1
+
+        u = g.u
+
+        # this is the MC limiter from van Leer (1977), as given in
+        # LeVeque (2002).  Note that this is slightly different than
+        # the expression from Colella (1990)
+
+        dc = g.scratch_array()
+        dl = g.scratch_array()
+        dr = g.scratch_array()
+
+        dc[ib:ie+1] = 0.5*(u[ib+1:ie+2] - u[ib-1:ie  ])
+        dl[ib:ie+1] = u[ib+1:ie+2] - u[ib  :ie+1]
+        dr[ib:ie+1] = u[ib  :ie+1] - u[ib-1:ie  ]
+
+        # these where's do a minmod()
+        d1 = 2.0*np.where(np.fabs(dl) < np.fabs(dr), dl, dr)
+        d2 = np.where(np.fabs(dc) < np.fabs(d1), dc, d1)
+        ldeltau = np.where(dl*dr > 0.0, d2, 0.0)
+
+        # now the interface states.  Note that there are 1 more interfaces
+        # than zones
+        ul = g.scratch_array()
+        ur = g.scratch_array()
+
+        # are these indices right?
+        #
+        #  --+-----------------+------------------+
+        #     ^       i       ^ ^        i+1
+        #     ur(i)     ul(i+1) ur(i+1)
+        #
+        ur[ib:ie+2] = u[ib:ie+2] - \
+                      0.5*(1.0 + u[ib:ie+2]*dt/self.grid.dx)*ldeltau[ib:ie+2]
+
+        ul[ib+1:ie+2] = u[ib:ie+1] + \
+                        0.5*(1.0 - u[ib:ie+1]*dt/self.grid.dx)*ldeltau[ib:ie+1]
+
+        return ul, ur
+
+
+    def riemann(self, ul, ur):
+        """
+        Riemann problem for Burgers' equation.
+        """
+
+        S = 0.5*(ul + ur)
+        ushock = np.where(S > 0.0, ul, ur)
+        ushock = np.where(S == 0.0, 0.0, ushock)
+
+        # rarefaction solution
+        urare = np.where(ur <= 0.0, ur, 0.0)
+        urare = np.where(ul >= 0.0, ul, urare)
+
+        us = np.where(ul > ur, ushock, urare)
+
+        return 0.5*us*us
+
+
+    def update(self, dt, flux):
+        """ conservative update """
+
+        g = self.grid
+
+        unew = g.scratch_array()
+
+        unew[g.ilo:g.ihi+1] = g.u[g.ilo:g.ihi+1] + \
+            dt/g.dx * (flux[g.ilo:g.ihi+1] - flux[g.ilo+1:g.ihi+2])
+
+        return unew
+
+
+    def evolve(self, C, tmax):
+
+        self.t = 0.0
+
+        g = self.grid
+
+        # main evolution loop
+        while (self.t < tmax):
+
+            # fill the boundary conditions
+            g.fill_BCs()
+
+            # get the timestep
+            dt = self.timestep(C)
+
+            if (self.t + dt > tmax):
+                dt = tmax - self.t
+
+            # get the interface states
+            ul, ur = self.states(dt)
+
+            # solve the Riemann problem at all interfaces
+            flux = self.riemann(ul, ur)
+
+            # do the conservative update
+            unew = self.update(dt, flux)
+
+            self.grid.u[:] = unew[:]
+
+            self.t += dt
+    
+
+if __name__ == "__main__":
+    #-----------------------------------------------------------------------------
+    # sine
+
+    xmin = 0.0
+    xmax = 1.0
+    nx = 256
+    ng = 2
+    g = Grid1d(nx, ng, bc="periodic")
+
+    # maximum evolution time based on period for unit velocity
+    tmax = (xmax - xmin)/1.0
+
+    C = 0.8
+
+    plt.clf()
+    
+    s = Simulation(g)
+
+
+    for i in range(0, 10):
+        tend = (i+1)*0.02*tmax
+        s.init_cond("sine")
+
+        uinit = s.grid.u.copy()
+
+        s.evolve(C, tend)
+
+        c = 1.0 - (0.1 + i*0.1)
+        g = s.grid
+        plt.plot(g.x[g.ilo:g.ihi+1], g.u[g.ilo:g.ihi+1], color=str(c))
+
+
+    #print(f"Completed Execution in {timeTaken} seconds")
+
+    g = s.grid
+    plt.plot(g.x[g.ilo:g.ihi+1], uinit[g.ilo:g.ihi+1], ls=":", color="0.9", zorder=-1)
+
+    plt.xlabel("$x$")
+    plt.ylabel("$u$")
+    plt.savefig("fv-burger-sine.pdf")
+
+
+    #-----------------------------------------------------------------------------
+    # rarefaction
+    
+    xmin = 0.0
+    xmax = 1.0
+    nx = 256
+    ng = 2
+    g = Grid1d(nx, ng, bc="outflow")
+
+    # maximum evolution time based on period for unit velocity
+    tmax = (xmax - xmin)/1.0
+
+    C = 0.8
+
+    plt.clf()
+
+    s = Simulation(g)
+
+    for i in range(0, 10):
+        tend = (i+1)*0.02*tmax
+
+        s.init_cond("rarefaction")
+
+        uinit = s.grid.u.copy()
+
+        s.evolve(C, tend)
+
+        c = 1.0 - (0.1 + i*0.1)
+        plt.plot(g.x[g.ilo:g.ihi+1], g.u[g.ilo:g.ihi+1], color=str(c))
+
+
+    plt.plot(g.x[g.ilo:g.ihi+1], uinit[g.ilo:g.ihi+1], ls=":", color="0.9", zorder=-1)
+
+    plt.xlabel("$x$")
+    plt.ylabel("$u$")
+
+    plt.savefig("fv-burger-rarefaction.pdf")
+    
diff --git a/Python Burgers example/python-us.txt b/Python Burgers example/python-us.txt
new file mode 100644
index 0000000..d6118fc
--- /dev/null
+++ b/Python Burgers example/python-us.txt	
@@ -0,0 +1,256 @@
+0.999999999999999
+0.9999999999999991
+0.9999999999999992
+0.9999999999999993
+0.9999999999999994
+0.9999999999999996
+0.9999999999999997
+0.9999999999999998
+0.9999999999999999
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0
+1.0000000000000002
+1.0000000000000024
+1.0000000000000235
+1.0000000000001896
+1.0000000000014284
+1.0000000000104186
+1.0000000000740106
+1.0000000005137009
+1.0000000034954541
+1.0000000234001618
+1.0000001547117818
+1.000001013611244
+1.0000065872150568
+1.0000424961094072
+1.0002748504319279
+1.001743354167995
+1.0065958193203346
+1.014635990126267
+1.0249968395795583
+1.0368531572695936
+1.049550849499212
+1.0626478439177824
+1.075873722895156
+1.0890728851076001
+1.1021642864183863
+1.1151150564243955
+1.1279229999395857
+1.1406042755763512
+1.1531836461867673
+1.165686557726796
+1.17813347740903
+1.1905369175727507
+1.2029010201645876
+1.2152230526504362
+1.2274958713795525
+1.2397104117165654
+1.251857539092869
+1.2639289952620607
+1.2759175266673788
+1.2878164849050542
+1.2996192296523668
+1.3113185874253426
+1.322906491864095
+1.3343738127847697
+1.3457103073318923
+1.356904601717775
+1.3679441164604864
+1.3788148574446542
+1.389501003093958
+1.3999842171881165
+1.4102426223451991
+1.420250147029987
+1.429979624854128
+1.4393984135314737
+1.448401240104463
+1.4571037924348336
+1.4668068032924166
+1.4726561228919206
+1.3706864195203843
+0.6973878274598113
+0.5267488145231821
+0.531733371064891
+0.541571102461252
+0.5503190861848413
+0.5592282938096137
+0.568593965716811
+0.5782851489961105
+0.5882572217462919
+0.5984859850833824
+0.6089441977509232
+0.6196083129738961
+0.6304591460041207
+0.6414804734405799
+0.6526580901517363
+0.6639792657296127
+0.6754325022401995
+0.6870075005015289
+0.6986951970381017
+0.7104877226393795
+0.722378149215317
+0.7343599509815423
+0.7464262255292186
+0.7585688858202705
+0.7707781922571594
+0.783043064764136
+0.795352520791625
+0.807698287777906
+0.8200781725764194
+0.8324992662618946
+0.8449797314428379
+0.8575479648727261
+0.8702384291711712
+0.883084257119224
+0.8961073767093505
+0.9093068160962303
+0.9226449261835163
+0.9360293298806136
+0.9492857489626352
+0.9621192294959894
+0.9740717054009133
+0.9844933811154147
+0.9926082003200034
+0.9977310106981597
+0.9995958908369349
+0.9999290616935412
+0.9999876088492843
+0.999997858268876
+0.9999996340859503
+0.9999999384058849
+0.999999989824941
+0.9999999983567549
+0.9999999997415411
+0.999999999960553
+0.9999999999941773
+0.9999999999991708
+0.9999999999998853
+0.9999999999999832
+0.9999999999999957
+0.9999999999999969
+0.9999999999999971
+0.9999999999999972
+0.9999999999999973
+0.9999999999999974
+0.9999999999999976
+0.9999999999999977
+0.9999999999999978
+0.9999999999999979
+0.999999999999998
+0.9999999999999981
+0.9999999999999982
+0.9999999999999983
+0.9999999999999984
+0.9999999999999986
+0.9999999999999987
+0.9999999999999988
+0.9999999999999989
diff --git a/Python Burgers example/python-xs.txt b/Python Burgers example/python-xs.txt
new file mode 100644
index 0000000..51b0e0d
--- /dev/null
+++ b/Python Burgers example/python-xs.txt	
@@ -0,0 +1,256 @@
+0.001953125
+0.005859375
+0.009765625
+0.013671875
+0.017578125
+0.021484375
+0.025390625
+0.029296875
+0.033203125
+0.037109375
+0.041015625
+0.044921875
+0.048828125
+0.052734375
+0.056640625
+0.060546875
+0.064453125
+0.068359375
+0.072265625
+0.076171875
+0.080078125
+0.083984375
+0.087890625
+0.091796875
+0.095703125
+0.099609375
+0.103515625
+0.107421875
+0.111328125
+0.115234375
+0.119140625
+0.123046875
+0.126953125
+0.130859375
+0.134765625
+0.138671875
+0.142578125
+0.146484375
+0.150390625
+0.154296875
+0.158203125
+0.162109375
+0.166015625
+0.169921875
+0.173828125
+0.177734375
+0.181640625
+0.185546875
+0.189453125
+0.193359375
+0.197265625
+0.201171875
+0.205078125
+0.208984375
+0.212890625
+0.216796875
+0.220703125
+0.224609375
+0.228515625
+0.232421875
+0.236328125
+0.240234375
+0.244140625
+0.248046875
+0.251953125
+0.255859375
+0.259765625
+0.263671875
+0.267578125
+0.271484375
+0.275390625
+0.279296875
+0.283203125
+0.287109375
+0.291015625
+0.294921875
+0.298828125
+0.302734375
+0.306640625
+0.310546875
+0.314453125
+0.318359375
+0.322265625
+0.326171875
+0.330078125
+0.333984375
+0.337890625
+0.341796875
+0.345703125
+0.349609375
+0.353515625
+0.357421875
+0.361328125
+0.365234375
+0.369140625
+0.373046875
+0.376953125
+0.380859375
+0.384765625
+0.388671875
+0.392578125
+0.396484375
+0.400390625
+0.404296875
+0.408203125
+0.412109375
+0.416015625
+0.419921875
+0.423828125
+0.427734375
+0.431640625
+0.435546875
+0.439453125
+0.443359375
+0.447265625
+0.451171875
+0.455078125
+0.458984375
+0.462890625
+0.466796875
+0.470703125
+0.474609375
+0.478515625
+0.482421875
+0.486328125
+0.490234375
+0.494140625
+0.498046875
+0.501953125
+0.505859375
+0.509765625
+0.513671875
+0.517578125
+0.521484375
+0.525390625
+0.529296875
+0.533203125
+0.537109375
+0.541015625
+0.544921875
+0.548828125
+0.552734375
+0.556640625
+0.560546875
+0.564453125
+0.568359375
+0.572265625
+0.576171875
+0.580078125
+0.583984375
+0.587890625
+0.591796875
+0.595703125
+0.599609375
+0.603515625
+0.607421875
+0.611328125
+0.615234375
+0.619140625
+0.623046875
+0.626953125
+0.630859375
+0.634765625
+0.638671875
+0.642578125
+0.646484375
+0.650390625
+0.654296875
+0.658203125
+0.662109375
+0.666015625
+0.669921875
+0.673828125
+0.677734375
+0.681640625
+0.685546875
+0.689453125
+0.693359375
+0.697265625
+0.701171875
+0.705078125
+0.708984375
+0.712890625
+0.716796875
+0.720703125
+0.724609375
+0.728515625
+0.732421875
+0.736328125
+0.740234375
+0.744140625
+0.748046875
+0.751953125
+0.755859375
+0.759765625
+0.763671875
+0.767578125
+0.771484375
+0.775390625
+0.779296875
+0.783203125
+0.787109375
+0.791015625
+0.794921875
+0.798828125
+0.802734375
+0.806640625
+0.810546875
+0.814453125
+0.818359375
+0.822265625
+0.826171875
+0.830078125
+0.833984375
+0.837890625
+0.841796875
+0.845703125
+0.849609375
+0.853515625
+0.857421875
+0.861328125
+0.865234375
+0.869140625
+0.873046875
+0.876953125
+0.880859375
+0.884765625
+0.888671875
+0.892578125
+0.896484375
+0.900390625
+0.904296875
+0.908203125
+0.912109375
+0.916015625
+0.919921875
+0.923828125
+0.927734375
+0.931640625
+0.935546875
+0.939453125
+0.943359375
+0.947265625
+0.951171875
+0.955078125
+0.958984375
+0.962890625
+0.966796875
+0.970703125
+0.974609375
+0.978515625
+0.982421875
+0.986328125
+0.990234375
+0.994140625
+0.998046875
diff --git a/Python Burgers example/run-burgers.py b/Python Burgers example/run-burgers.py
new file mode 100644
index 0000000..a4d1dbc
--- /dev/null
+++ b/Python Burgers example/run-burgers.py	
@@ -0,0 +1,20 @@
+import benchmarkedBurgers as burgers
+import math
+
+totalTime = 0
+nTests = 3
+results = []
+for i in range(nTests):
+    time = burgers.main()
+    results.append(time)
+    totalTime += time
+avg = totalTime / nTests
+totaldif = 0
+for i in results:
+    totaldif += (i - avg)**2
+sd = math.sqrt(totaldif / (totalTime - 1))
+
+avgms = avg * 1000
+print("Average time: " + str(avg))
+print("In Milliseconds: " + str(avgms))
+print("SD: " + str(sd))
diff --git a/collectedGraphs/Fluxes/TrixiFluxEngquistOsher.png b/collectedGraphs/Fluxes/TrixiFluxEngquistOsher.png
new file mode 100644
index 0000000..5058168
Binary files /dev/null and b/collectedGraphs/Fluxes/TrixiFluxEngquistOsher.png differ
diff --git a/collectedGraphs/Fluxes/TrixiFluxEngquistOsherComparison.png b/collectedGraphs/Fluxes/TrixiFluxEngquistOsherComparison.png
new file mode 100644
index 0000000..6add732
Binary files /dev/null and b/collectedGraphs/Fluxes/TrixiFluxEngquistOsherComparison.png differ
diff --git a/collectedGraphs/Fluxes/TrixiFluxGodunov.png b/collectedGraphs/Fluxes/TrixiFluxGodunov.png
new file mode 100644
index 0000000..a7abe52
Binary files /dev/null and b/collectedGraphs/Fluxes/TrixiFluxGodunov.png differ
diff --git a/collectedGraphs/Fluxes/TrixiFluxGodunovComparison.png b/collectedGraphs/Fluxes/TrixiFluxGodunovComparison.png
new file mode 100644
index 0000000..63088e3
Binary files /dev/null and b/collectedGraphs/Fluxes/TrixiFluxGodunovComparison.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly2.png b/collectedGraphs/PolynomialDegree/TrixiPoly2.png
new file mode 100644
index 0000000..d32955f
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly2.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly2Comparison.png b/collectedGraphs/PolynomialDegree/TrixiPoly2Comparison.png
new file mode 100644
index 0000000..c06280d
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly2Comparison.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly3.png b/collectedGraphs/PolynomialDegree/TrixiPoly3.png
new file mode 100644
index 0000000..0583bec
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly3.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly3Comparison.png b/collectedGraphs/PolynomialDegree/TrixiPoly3Comparison.png
new file mode 100644
index 0000000..0583bec
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly3Comparison.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly4.png b/collectedGraphs/PolynomialDegree/TrixiPoly4.png
new file mode 100644
index 0000000..eaae35e
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly4.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly4Comp.png b/collectedGraphs/PolynomialDegree/TrixiPoly4Comp.png
new file mode 100644
index 0000000..26a1ee1
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly4Comp.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly5.png b/collectedGraphs/PolynomialDegree/TrixiPoly5.png
new file mode 100644
index 0000000..c015b02
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly5.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly5Comp.png b/collectedGraphs/PolynomialDegree/TrixiPoly5Comp.png
new file mode 100644
index 0000000..1bf12b8
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly5Comp.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly6.png b/collectedGraphs/PolynomialDegree/TrixiPoly6.png
new file mode 100644
index 0000000..b300f08
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly6.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly6Comp.png b/collectedGraphs/PolynomialDegree/TrixiPoly6Comp.png
new file mode 100644
index 0000000..f9d0c93
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly6Comp.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly7.png b/collectedGraphs/PolynomialDegree/TrixiPoly7.png
new file mode 100644
index 0000000..71d388c
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly7.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly7Comp.png b/collectedGraphs/PolynomialDegree/TrixiPoly7Comp.png
new file mode 100644
index 0000000..4630d5d
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly7Comp.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly8.png b/collectedGraphs/PolynomialDegree/TrixiPoly8.png
new file mode 100644
index 0000000..4ffbb54
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly8.png differ
diff --git a/collectedGraphs/PolynomialDegree/TrixiPoly8Comp.png b/collectedGraphs/PolynomialDegree/TrixiPoly8Comp.png
new file mode 100644
index 0000000..37e9fb6
Binary files /dev/null and b/collectedGraphs/PolynomialDegree/TrixiPoly8Comp.png differ
diff --git a/collectedGraphs/PythonVsTrixi/PythonVsTrixi-1024.png b/collectedGraphs/PythonVsTrixi/PythonVsTrixi-1024.png
new file mode 100644
index 0000000..fa64aa6
Binary files /dev/null and b/collectedGraphs/PythonVsTrixi/PythonVsTrixi-1024.png differ
diff --git a/collectedGraphs/PythonVsTrixi/PythonVsTrixi-256.png b/collectedGraphs/PythonVsTrixi/PythonVsTrixi-256.png
new file mode 100644
index 0000000..69bc92d
Binary files /dev/null and b/collectedGraphs/PythonVsTrixi/PythonVsTrixi-256.png differ
diff --git a/collectedGraphs/errors.txt b/collectedGraphs/errors.txt
new file mode 100644
index 0000000..2cbb62f
--- /dev/null
+++ b/collectedGraphs/errors.txt
@@ -0,0 +1,317 @@
+python 256 vs trixi 256 refinement 6
+dy: 0.0048056807547305624   
+abs_dy: 0.007653058515335133
+relerr: 0.008469341114883688
+pererr: 0.846934111488369
+mean_err: 0.007653058515335133
+MSE: 0.0018368798814202886
+
+p v t 1024 r 8 
+dy: -0.00039388333544266034
+abs_dy: 0.0011925272442834382
+relerr: 0.001191920299046318
+pererr: 0.11919202990463126
+mean_err: 0.0011925272442834382
+MSE: 0.0001374401231385753
+
+trixi init refinement tests - this doesn't work bc the values its testing are of course the ones that are correct
+1
+9138 samples
+238μs to 2.7s
+median 257μs
+mean 569μs
+43.90KiB
+
+dy: -0.06287360535080128
+abs_dy: 0.1405070061025938
+relerr: 0.17204119586133576
+pererr: 17.204119586133576
+mean_err: 0.1405070061025938
+MSE: 0.03731663559015569
+
+2
+median 381μs
+mean 770μs
+55.95KiB
+
+dy: 0.0782067306835732
+abs_dy: 0.09439203823313455
+relerr: 0.09174725477864251
+pererr: 9.174725477864254
+mean_err: 0.09439203823313455
+MSE: 0.027238624870018876
+
+3
+median 604μs
+1.355ms
+87.82KiB
+
+dy: 0.017511428188698958
+abs_dy: 0.03617746939321864
+relerr: 0.04552322470740188
+pererr: 4.552322470740188
+mean_err: 0.03617746939321864
+MSE: 0.00557124123320567
+
+4
+2758 samples
+median 928μs
+2.3301ms
+163.27KiB
+
+dy: -0.019213917632710673
+abs_dy: 0.03370475679173245
+relerr: 0.04254520835315101
+pererr: 4.254520835315101
+mean_err: 0.03370475679173245
+MSE: 0.01882122012780893
+
+5
+1124 samples
+median 1.959ms
+mean 5.210ms
+411.91KiB
+
+dy: -0.006428145548389685
+abs_dy: 0.013606207256852229
+relerr: 0.015855393992331483
+pererr: 1.585539399233149
+mean_err: 0.013606207256852229
+MSE: 0.0036706011005593073
+
+6
+943 samples
+4.879ms median
+mean 5.295ms
+1.19MiB
+
+dy: 0.0048056807547305624
+abs_dy: 0.007653058515335133
+relerr: 0.008469341114883688
+pererr: 0.846934111488369
+mean_err: 0.007653058515335133
+MSE: 0.0018368798814202886
+
+7
+357 samples
+median 13.692ms
+mean 14ms
+3.82MiB
+
+256
+dy: -0.0003859279601094814
+abs_dy: 0.0034390741592044637
+relerr: 0.003984757033531754
+pererr: 0.3984757033531749
+mean_err: 0.0034390741592044637
+MSE: 0.0006419785617369417
+
+8
+116 samples
+median 42ms
+mean 43.32ms
+13.02MiB
+
+256
+dy: -0.0012501483748248495
+abs_dy: 0.002576959312279194
+relerr: 0.0029673821928554714
+pererr: 0.29673821928554733
+mean_err: 0.002576959312279194
+MSE: 0.00010065629134159438
+
+1024
+dy: -0.00039388333544266034
+abs_dy: 0.0011925272442834382
+relerr: 0.001191920299046318
+pererr: 0.11919202990463126
+mean_err: 0.0011925272442834382
+MSE: 0.000137440123138575
+
+9
+22 samples
+median 124ms
+mean 277ms
+44.44 MiB
+
+256
+dy: -0.0017425669460848104
+abs_dy: 0.0031270411105882414
+relerr: 0.003442321994161419
+pererr: 0.34423219941614236
+mean_err: 0.0031270411105882414
+MSE: 0.00021399684952601075
+
+10
+6 samples
+median 627ms
+mean 849ms
+161.08MiB
+
+256
+dy: -0.0020377351469643462
+abs_dy: 0.0034214312170370823
+relerr: 0.003657786243589712
+pererr: 0.36577862435897146
+mean_err: 0.0034214312170370823
+MSE: 0.0003252164958825785
+
+
+**Trixi vs Julia**
+-Julia simple method
+
+-Julia spectral
+
+-Trixi
+
+ 
+N 1024 for all
+~Surface flux/DGSEM method~
+-flux lax friedrichs
+:what we have been using so far so refer to prev results
+
+-flux godunov
+148 samples
+33.791 ms (+- 6.178ms (1 SD))
+11.76 MiB
+
+dy: -0.00037428645564645485
+abs_dy: 0.0010918283812354091
+relerr: 0.001114385434592326
+pererr: 0.11143854345923196
+mean_err: 0.0010918283812354091
+MSE: 0.0001220728747760191
+
+-flux engquist osher
+146 samples
+34.287 ms
+11.76 MiB
+
+dy: -0.00037428645564645485
+abs_dy: 0.0010918283812354091
+relerr: 0.001114385434592326
+pererr: 0.11143854345923196
+mean_err: 0.0010918283812354091
+MSE: 0.0001220728747760191
+
+(these are the ones allowed by Trixi in Burgers, found by looking at source code)
+-flux_ranocha
+error
+
+used 1024
+~Polynomial degree~
+2
+range 16.8 ms to 51.3 ms
+median 18 ms
+mean 18.7 ms
+6.35 MiB
+
+dy: 4.472588937928038e-5
+abs_dy: 0.0014646496510331952
+relerr: 0.0013328690228656482
+pererr: 0.13328690228656417
+mean_err: 0.0014646496510331952
+MSE: 0.0001279568437868768
+
+3
+see other results
+
+4
+55ms to 82 ms
+median 57.7
+mean 60.35ms
+21.92 MiB
+
+dy: -0.0013089507679571355
+abs_dy: 0.0029887766189824632
+relerr: 0.003431487446766604
+pererr: 0.34314874467666034
+mean_err: 0.0029887766189824632
+MSE: 0.00020116362338661384
+
+5
+86ms to 115ms
+median 89.8ms
+mean 93.53ms
+33.20MiB
+
+dy: -0.0005780374942881585
+abs_dy: 0.0008673755379447632
+relerr: 0.0009720248845475113
+pererr: 0.09720248845475078
+mean_err: 0.0008673755379447632
+MSE: 0.0003736103229096141
+
+6
+136ms to 179ms
+median 153ms
+man 152.26ms
+49.67 MiB
+
+dy: -0.0006139748444910812
+abs_dy: 0.0010852029188884898
+relerr: 0.0012005977225163333
+pererr: 0.12005977225163275
+mean_err: 0.0010852029188884898
+MSE: 0.0006104670474898316
+
+7
+201ms to 1.8s
+226ms median
+426ms mean
+69.11MiB
+
+dy: -0.00022281978255088701
+abs_dy: 0.0005272251879351443
+relerr: 0.0004942344931976489
+pererr: 0.0494234493197646
+mean_err: 0.0005272251879351443
+MSE: 8.717804097281288e-5
+
+8
+9 samples
+320ms to 2s
+median 362ms
+mean 590ms
+
+dy: -0.0002972028112626587
+abs_dy: 0.00048692105057538207
+relerr: 0.0004697981984329635
+pererr: 0.046979819843296085
+mean_err: 0.00048692105057538207
+MSE: 0.00011014404909550833
+
+89.91MiB
+
+
+
+results
+python (N = 256) (127 samples) - 53.03 ms
+python (N = 1024) (127 samples) - 271.8 ms
+julia (N = 256, refinement 4) (3256 samples) - 2.02 ms, 175 kb
+julia (N = 1024, init refinement 8) (128 samples) - 39.07 ms, 13 mb
+julia (N = 1024, init refinement 9) (28 samples) - 201.47 ms
+julia (N = 16384, init refinement 10) (10 samples) - 537.77 ms
+
+NEW PYTHON RESULTS (not including plotting time)
+python (N = 256) (127 samples) - 37.66 ms
+python (N = 1024) (127 samples) - 260.22 ms
+python (N = 256) (3256 samples) - 40.34 ms
+
+NEW JULIA RESULTS
+julia (N = 256, refinement 5) (1436 samples) - 4.83 ms, 420 kb
+julia (N = 256, refinement 6) (1033 samples) - 4.83 ms, 1.19 mb
+refinement 6 is true 256!!
+julia (N = 256, refinement 7) (368 samples) - 13.57 ms, 3.81 mb
+
+WITH NO PLOTTING
+n = 8 1.18ms
+n 16 2.03ms
+n 32 3.92ms
+n 64 7.60ms
+n 128 17ms
+n 256 37.79ms
+n 512 95.25ms
+1024 254ms
+n 2048 797 ms
+n 4096 2.56s (2562ms)
\ No newline at end of file
diff --git a/collectedGraphs/initRefinementTests/InitRefinement1.png b/collectedGraphs/initRefinementTests/InitRefinement1.png
new file mode 100644
index 0000000..156e8f0
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement1.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement10Compare.png b/collectedGraphs/initRefinementTests/InitRefinement10Compare.png
new file mode 100644
index 0000000..113a215
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement10Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement1Compare.png b/collectedGraphs/initRefinementTests/InitRefinement1Compare.png
new file mode 100644
index 0000000..e825000
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement1Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement2.png b/collectedGraphs/initRefinementTests/InitRefinement2.png
new file mode 100644
index 0000000..677eccb
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement2.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement2Compare.png b/collectedGraphs/initRefinementTests/InitRefinement2Compare.png
new file mode 100644
index 0000000..1695b5f
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement2Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement3.png b/collectedGraphs/initRefinementTests/InitRefinement3.png
new file mode 100644
index 0000000..57632a6
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement3.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement3Compare.png b/collectedGraphs/initRefinementTests/InitRefinement3Compare.png
new file mode 100644
index 0000000..291a6d3
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement3Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement4.png b/collectedGraphs/initRefinementTests/InitRefinement4.png
new file mode 100644
index 0000000..17d654f
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement4.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement4Compare.png b/collectedGraphs/initRefinementTests/InitRefinement4Compare.png
new file mode 100644
index 0000000..cc7be37
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement4Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement5.png b/collectedGraphs/initRefinementTests/InitRefinement5.png
new file mode 100644
index 0000000..8327ea6
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement5.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement5Compare.png b/collectedGraphs/initRefinementTests/InitRefinement5Compare.png
new file mode 100644
index 0000000..d866f38
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement5Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement6.png b/collectedGraphs/initRefinementTests/InitRefinement6.png
new file mode 100644
index 0000000..386e9b3
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement6.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement6Compare.png b/collectedGraphs/initRefinementTests/InitRefinement6Compare.png
new file mode 100644
index 0000000..7b5f037
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement6Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement7.png b/collectedGraphs/initRefinementTests/InitRefinement7.png
new file mode 100644
index 0000000..0c29e84
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement7.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement7Compare.png b/collectedGraphs/initRefinementTests/InitRefinement7Compare.png
new file mode 100644
index 0000000..9d2b173
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement7Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement8.png b/collectedGraphs/initRefinementTests/InitRefinement8.png
new file mode 100644
index 0000000..4a9012a
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement8.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement8Compare.png b/collectedGraphs/initRefinementTests/InitRefinement8Compare.png
new file mode 100644
index 0000000..52c2034
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement8Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement9.png b/collectedGraphs/initRefinementTests/InitRefinement9.png
new file mode 100644
index 0000000..0f36201
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement9.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinement9Compare.png b/collectedGraphs/initRefinementTests/InitRefinement9Compare.png
new file mode 100644
index 0000000..24a50b1
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinement9Compare.png differ
diff --git a/collectedGraphs/initRefinementTests/InitRefinementsOverlaid.png b/collectedGraphs/initRefinementTests/InitRefinementsOverlaid.png
new file mode 100644
index 0000000..2ccda8a
Binary files /dev/null and b/collectedGraphs/initRefinementTests/InitRefinementsOverlaid.png differ
diff --git a/collectedGraphs/initRefinementTests/r1.gif b/collectedGraphs/initRefinementTests/r1.gif
new file mode 100644
index 0000000..cd42016
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r1.gif differ
diff --git a/collectedGraphs/initRefinementTests/r2.gif b/collectedGraphs/initRefinementTests/r2.gif
new file mode 100644
index 0000000..35c73c3
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r2.gif differ
diff --git a/collectedGraphs/initRefinementTests/r3.gif b/collectedGraphs/initRefinementTests/r3.gif
new file mode 100644
index 0000000..05ba578
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r3.gif differ
diff --git a/collectedGraphs/initRefinementTests/r4.gif b/collectedGraphs/initRefinementTests/r4.gif
new file mode 100644
index 0000000..c6cdfbc
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r4.gif differ
diff --git a/collectedGraphs/initRefinementTests/r5.gif b/collectedGraphs/initRefinementTests/r5.gif
new file mode 100644
index 0000000..6300df3
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r5.gif differ
diff --git a/collectedGraphs/initRefinementTests/r6.gif b/collectedGraphs/initRefinementTests/r6.gif
new file mode 100644
index 0000000..9c878b1
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r6.gif differ
diff --git a/collectedGraphs/initRefinementTests/r7.gif b/collectedGraphs/initRefinementTests/r7.gif
new file mode 100644
index 0000000..06d275d
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r7.gif differ
diff --git a/collectedGraphs/initRefinementTests/r8.gif b/collectedGraphs/initRefinementTests/r8.gif
new file mode 100644
index 0000000..2870bc6
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r8.gif differ
diff --git a/collectedGraphs/initRefinementTests/r9.gif b/collectedGraphs/initRefinementTests/r9.gif
new file mode 100644
index 0000000..4ce15f3
Binary files /dev/null and b/collectedGraphs/initRefinementTests/r9.gif differ
diff --git a/examples/burger-chef.jl b/examples/burger-chef.jl
index 0e92cb3..44eff07 100644
--- a/examples/burger-chef.jl
+++ b/examples/burger-chef.jl
@@ -4,26 +4,56 @@ using DiscSimulations
 using Plots
 using DifferentialEquations
 using Printf
+using Trixi
 
 N = 512
 xmax = 2π
 
-simple_x, _, simple_problem = BurgerSimple.setup(N, xmax)
+function initial_condition_burgers(x)
+    u = 1.0
+    if(x[1] >= 0.333 && x[1] <= 0.666)
+        u = u + 0.5*sin(2.0*π*(x[1]-0.333)/0.333)
+    end
+    return u
+end
+
+function _wrap_init_conditions(f::Function)
+    function _wrapper(x)
+        map(x) do i
+            f(SVector(i))
+        end
+    end
+end
+
+simple_x, _, simple_problem = BurgerSimple.setup(N, xmax, _wrap_init_conditions(initial_condition_burgers))
 simple_sol =
     @time solve(simple_problem, Rodas5(autodiff = false); reltol = 1e-7, abstol = 1e-7)
 
-spectral_x, _, T, Ti, spectral_problem = BurgerSpectral.setup(N, xmax)
+spectral_x, _, T, Ti, spectral_problem = BurgerSpectral.setup(N, xmax, _wrap_init_conditions(initial_condition_burgers))
 spectral_sol =
     @time solve(spectral_problem, Rodas5(autodiff = false); reltol = 1e-7, abstol = 1e-7)
 
-# trixi_x, trixi_problem = BurgerTrixi.setup(N, xmax)
-# trixi_sol = @time solve(trixi_problem, RDPK3SpFSAL49(), abstol=1.0e-7, reltol=1.0e-7)
+trixi_x, trixi_problem = BurgerTrixi.setup(0, xmax, initial_condition_burgers, (0.0,2.0), 7, 3, flux_lax_friedrichs)
+trixi_sol = @time solve(trixi_problem, RDPK3SpFSAL49(), abstol=1.0e-7, reltol=1.0e-7)
 
 # use the new disc solution type
-discsol = BurgerTrixi.solve_disc(N, xmax)
+discsol = BurgerTrixi.solve_disc(0, xmax, 7)
 # plot that individually
 DiscSimulations.plotgif(discsol, 0.0, 1.0)
 
+p = plot()
+plot!(simple_x, simple_sol(1.0), label = "Simple")
+plot!(spectral_x, Ti * spectral_sol(1.0), label = "Spectral")
+plot!(trixi_x, Ti * trixi_sol(1.0), label = "Trixi")
+plot!(
+    p,
+    legend = :outertopright,
+    title = Printf.@sprintf("t = %1.2f", 1.0),
+    ylims = (-0.1, 2.1),
+    xlabel = "x",
+    ylabel = "f",
+)
+
 tspan = collect(range(0.0, 1.0, 60))
 frames = Plots.@animate for t in tspan
     p = plot()
diff --git a/examples/younsi-2012.jl b/examples/younsi-2012.jl
deleted file mode 100644
index 40fc411..0000000
--- a/examples/younsi-2012.jl
+++ /dev/null
@@ -1,25 +0,0 @@
-using DiscSimulations
-using Plots
-
-# black hole mass
-M = 1.0
-
-# black hole spin
-a = 0.2
-
-# keplerian radius
-rₖ = 12.0
-
-# angular velocity profile index of the torus
-n = 0.21
-
-begin
-    p = plot(xlabel = "r", ylabel = "z", legend = :outertopright)
-    for a in [0.0, 0.2, 0.4, 0.6, 0.8]
-        r, z = @time isobar(M, a, n, rₖ)
-        plot!(r, z, label = "a = $a")
-        # plot but upside down
-        plot!(r, -z, label = "a = $a")
-    end
-    p
-end
diff --git a/src/BurgersTestCase.jl b/src/BurgersTestCase.jl
new file mode 100644
index 0000000..83e341d
--- /dev/null
+++ b/src/BurgersTestCase.jl
@@ -0,0 +1,112 @@
+include("DiscSimulations.jl")
+using .DiscSimulations
+using BenchmarkTools
+using Trixi
+using Plots
+
+xmin = 0.0
+xmax = 1.0
+
+tmax = (xmax - xmin)/1.0
+t_span = (0.0, 0.2*tmax)
+
+function initial_condition_burgers(x)
+    u = 1.0
+    if(x[1] >= 0.333 && x[1] <= 0.666)
+        u = u + 0.5*sin(2.0*π*(x[1]-0.333)/0.333)
+    end
+    return u
+end
+
+#params expects this so we're defining a dummy one
+source_zero(u, x, t, equations) = SVector(0, 0, 0)
+
+params = DiscSimulations.Parameters(xmin, xmax, initial_condition_burgers, source_zero, t_span, 8, 3, Trixi.flux_lax_friedrichs)
+
+@benchmark solution = DiscSimulations.main(params, DiscSimulations.BurgersSimulation())
+solution = DiscSimulations.main(params, DiscSimulations.BurgersSimulation())
+s = DiscSimulations.DiscSolution(solution.sol, solution.semi, DiscSimulations.OneDimension())
+
+function importAndRead(filexs, fileus)
+    pyDXs = []
+    pyDUs = []
+ 
+    #import python data
+    pyDataX = open(filexs, "r")
+    line = 0  
+    while ! eof(pyDataX)          
+       s = parse(Float64, readline(pyDataX))         
+       line += 1
+       push!(pyDXs, s)
+    end #while
+    close(pyDataX)
+ 
+    pyDataU = open(fileus, "r")
+    line = 0  
+    while ! eof(pyDataU)          
+       s = parse(Float64, readline(pyDataU)) 
+       line += 1
+       push!(pyDUs, s)
+    end #while
+    close(pyDataU)
+ 
+    return pyDXs, pyDUs
+ end 
+
+pyDXs1024, pyDUs1024 = importAndRead("python-xs-1024.txt", "python-us-1024.txt")
+pyDXs256, pyDUs256 = importAndRead("python-xs-256.txt", "python-us-256.txt")
+
+pd = PlotData1D(solution.sol)
+
+plot(pd, title="Trixi and Python Burgers Comparison", ylabel="u")
+plot!(pyDXs, pyDUs)
+savefig("PythonVsTrixi-1024.png")
+
+#error calculation
+t_dy = 0
+t_abs_dy = 0   
+t_relerr = 0  
+t_pererr = 0   
+t_mean_err = 0    
+t_MSE = 0        
+t_RMSE = 0
+
+y = last(poly4.sol)
+inc = length(y) / 256
+inc = 256 / length(y)
+for i in 1:256
+    #y0 = pyDUs256[trunc(Int, i*inc)]
+    #y1 = y[i]
+
+    y0 = pyDUs256[i]
+    y1 = y[trunc(Int, i*inc)]
+
+    dy = y0-y1 # error 
+    abs_dy = abs(y0-y1)   # absolute error 
+    relerr = abs(y0-y1)./y0  # relative error 
+    pererr = abs(y0-y1)./y0*100   # percentage error 
+    mean_err = mean(abs(y0-y1))    # mean absolute error 
+    MSE = mean((y0-y1).^2)        # Mean square error 
+    #RMSE = sqrt(mean((y0-y1).^2)) # Root mean square error
+
+    t_dy += dy
+    t_abs_dy += abs_dy   
+    t_relerr += relerr  
+    t_pererr += pererr   
+    t_mean_err += mean_err    
+    t_MSE += MSE        
+    #t_RMSE += RSME
+end
+
+function printErrs(t_dy, t_abs_dy, t_relerr, t_pererr, t_mean_err, t_MSE, N)
+    print("dy: ", string(t_dy / N), "\n")
+    print("abs_dy: ", string(t_abs_dy / N), "\n")   
+    print("relerr: ", string(t_relerr / N), "\n")  
+    print("pererr: ", string(t_pererr / N), "\n")   
+    print("mean_err: ", string(t_mean_err / N), "\n")    
+    print("MSE: ", string(t_MSE / N), "\n")        
+end
+
+printErrs(t_dy, t_abs_dy, t_relerr, t_pererr, t_mean_err, t_MSE, 256)
+
+DiscSimulations.plotgif(poly2, t_span[1], t_span[2]; xlims = (0, 1), ylims = (0, 2), legend = :topright)
\ No newline at end of file
diff --git a/src/DiscSimulations.jl b/src/DiscSimulations.jl
index e113c97..da9f684 100644
--- a/src/DiscSimulations.jl
+++ b/src/DiscSimulations.jl
@@ -1,38 +1,60 @@
-module DiscSimulationsjl
+module DiscSimulations
 
 using LinearAlgebra
 using Plots
 using Printf
-
-include("solution.jl")
-
-include("burgers/utils.jl")
-#include("burgers/spectral.jl")
-#include("burgers/simple.jl")
-include("simpleGravitySim/simulation.jl")
-include("younsi-2012.jl")
-include("Euler_with_source.jl")
-include("burgers/trixi.jl")
-
-struct Parameters{}
-    N
-    xmin
-    xmax
+using TruncatedStacktraces
+using DifferentialEquations
+TruncatedStacktraces.VERBOSE[] = true
+
+#struct Parameters{T,DensityFuncType,SourceFuncType}
+struct Parameters{T}
+    xmin::T
+    xmax::T
     ρ
     source
+    t_span
+    refinement_level::Int
+    polydegree
+    flux_type
 end
 
 abstract type SimulationType end
-main(p::Parameters, ::SimulationType) = error("Unknown simulation type.")
+main(p::Parameters, ::SimulationType) = error("Unknown simulation type: $(typeof(t))")
 
 struct BurgersSimulation <: SimulationType end
 struct EulerWithSourceSimulation <: SimulationType end
 struct SimpleGravitySimulation <: SimulationType end
+struct SimpleBurgersSimulation <: SimulationType end
+struct SpectralBurgersSimulation <: SimulationType end
 
-main(p::Parameters, ::BurgersSimulation) = BurgerTrixi.solve_disc(p.N, p.xmax)
+main(p::Parameters, ::BurgersSimulation) = BurgerTrixi.solve_disc(p.xmin, p.xmax, p.refinement_level, p.ρ, p.t_span, p.polydegree, p.flux_type)
 main(p::Parameters, ::EulerWithSourceSimulation) = EulerSource.solve_euler(p.source)
 main(p::Parameters, ::SimpleGravitySimulation) = SimpleGravitySimulation.solve_gravity_sim(p.ρ)
 
+function main(p::Parameters, ::SimpleBurgersSimulation)
+    simple_x, _, simple_problem = BurgerSimple.setup(2^(p.refinement_level + 2), p.xmax, p.ρ, p.t_span)
+    simple_sol = @time solve(simple_problem, Rodas5(autodiff = false); reltol = 1e-7, abstol = 1e-7)
+    return simple_sol
+end
+function main(p::Parameters, ::SpectralBurgersSimulation)
+    spectral_x, _, T, Ti, spectral_problem = BurgerSpectral.setup(2^(p.refinement_level + 2), p.xmax, p.ρ, p.t_span)
+    spectral_sol = @time solve(spectral_problem, Rodas5(autodiff = false); reltol = 1e-7, abstol = 1e-7)
+    return spectral_sol
+end
+
+include("solution.jl")
+
+include("burgers/utils.jl")
+include("burgers/trixi.jl")
+include("burgers/spectral.jl")
+include("burgers/simple.jl")
+include("simpleGravitySim/simpleGravity.jl")
+include("simpleGravitySim/simulation.jl")
+#include("younsi-2012.jl")
+include("Euler_with_source.jl")
+
+
 export Parameters, SimulationType, BurgersSimulation, EulerWithSourceSimulation, SimpleGravitySimulation
 
 end # module DiscSimulations
diff --git a/src/Euler_with_source.jl b/src/Euler_with_source.jl
index 7382557..80c97d7 100644
--- a/src/Euler_with_source.jl
+++ b/src/Euler_with_source.jl
@@ -20,7 +20,7 @@ using Plots
     du1 = -cos(ϕ)*ρ*(G*mₛ/(r^2))
     du2 = -sin(ϕ)*ρ*(G*mₛ/(r^2))
     du3 = -(ux*cos(ϕ) + uy*sin(ϕ))*ρ*(G*mₛ/(r^2))
-    return SVector(du1, du2, du3)
+    return SVector(0, du1, du3)
 end
 
 source_zero(u, x, t, equations) = SVector(0, 0, 0)
diff --git a/src/burgers/trixi.jl b/src/burgers/trixi.jl
index d035a4c..f3392ae 100644
--- a/src/burgers/trixi.jl
+++ b/src/burgers/trixi.jl
@@ -5,32 +5,26 @@ using DifferentialEquations
 using Plots
 using Printf
 
-include("utils.jl")
-using .BurgerUtils
-
-include("../solution.jl")
-
-#import DiscSimulations:
-#    BurgerUtils.STANDARD_BURGER_INIT, BurgerUtils.STANDARD_BURGER_TSPAN, DiscSolution, OneDimension
-
-function setup(N, x_max, init = BurgerUtils.STANDARD_BURGER_INIT, t_span = BurgerUtils.STANDARD_BURGER_TSPAN)
-    x_min = 0.0
+import DiscSimulations:
+    STANDARD_BURGER_INIT, STANDARD_BURGER_TSPAN, DiscSolution, OneDimension
 
+function setup(x_min, x_max, init, t_span, refinement_level, polydegree, flux_type)
+    N = 2^(refinement_level + 2)
     equations = InviscidBurgersEquation1D()
-    solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+    solver = DGSEM(polydeg = polydegree, surface_flux = flux_type)
 
-    mesh = TreeMesh(x_min, x_max, initial_refinement_level = 3, n_cells_max = N)
+    mesh = TreeMesh(x_min, x_max, initial_refinement_level = refinement_level, n_cells_max = N)
     semi = SemidiscretizationHyperbolic(mesh, equations, (x, t, eqs) -> init(x), solver)
 
-    problem = semidiscretize(semi, t_span)
+    ode = semidiscretize(semi, t_span)
 
-    semi, problem
+    semi, ode
 end
 
-function solve_disc(N, x_max, init = BurgerUtils.STANDARD_BURGER_INIT, t_span = BurgerUtils.STANDARD_BURGER_TSPAN)
-    semi, problem = setup(N, x_max, init, t_span)
+function solve_disc(x_min, x_max, refinement_level, init = STANDARD_BURGER_INIT, t_span = STANDARD_BURGER_TSPAN, polydegree = 3, flux_type = flux_lax_friedrichs)
+    semi, problem = setup(x_min, x_max, init, t_span, refinement_level, polydegree, flux_type)
     sol = solve(problem, RDPK3SpFSAL49(), abstol = 1.0e-7, reltol = 1.0e-7)
-    return DiscSolution(sol, semi)
+    return DiscSolution(sol, semi, OneDimension())
 end
 
 end # module BurgerTrixi
diff --git a/src/burgers/utils.jl b/src/burgers/utils.jl
index b537f3f..eb8c32b 100644
--- a/src/burgers/utils.jl
+++ b/src/burgers/utils.jl
@@ -1,5 +1,3 @@
-module BurgerUtils
-
 function central_difference(Δu, N)
     lower = -ones(Float64, N - 1)
     mid = zeros(Float64, N)
@@ -11,8 +9,4 @@ function STANDARD_BURGER_INIT(x)
     @. 1 - cos(x)
 end
 
-const STANDARD_BURGER_TSPAN = (0.0, 1.0)
-
-end # module BurgerUtils
-
-export BurgerUtils
\ No newline at end of file
+const STANDARD_BURGER_TSPAN = (0.0, 1.0)
\ No newline at end of file
diff --git a/src/simpleGravitySim/simulation.jl b/src/simpleGravitySim/simulation.jl
index d32ff4d..ab0d228 100644
--- a/src/simpleGravitySim/simulation.jl
+++ b/src/simpleGravitySim/simulation.jl
@@ -1,7 +1,7 @@
 module simpleGravitySimulation
 
-include("./simpleGravity.jl")
-using .simpleGravity
+import DiscSimulations:
+    simpleGravity
 using Trixi
 using OrdinaryDiffEq
 using Plots
diff --git a/src/solution.jl b/src/solution.jl
index 64c9b7f..b942179 100644
--- a/src/solution.jl
+++ b/src/solution.jl
@@ -7,7 +7,6 @@ struct OneDimension <: AbstractDimension end
 struct DiscSolution{Dimension,SolutionType,DiscretiziationType}
     sol::SolutionType
     semi::DiscretiziationType
-    #to add: keyword args
 end
 
 function DiscSolution(sol, semi, dimension_type::AbstractDimension)
@@ -15,25 +14,27 @@ function DiscSolution(sol, semi, dimension_type::AbstractDimension)
 end
 
 #1D version of plot
-function dataplot(solution::DiscSolution{OneDimension})
-    return PlotData1D(solution.sol(t), solution.semi)
+function dataplot(solution::DiscSolution{<:OneDimension}, solver, t)
+    return PlotData1D(solver(t), solution.semi)
 end
 
 #2D version of plot
-function dataplot(solution::DiscSolution{TwoDimension})
-    return PlotData2D(solution.sol(t), solution.semi)
+function dataplot(solution::DiscSolution{<:TwoDimension}, solver, t)
+    return PlotData2D(solver(t), solution.semi)
 end
 
-function plotgif(solution, tmin, tmax)
+dataplot(solution::DiscSolution{<:AbstractDimension}, t) = error("Unknown simulation type: $(typeof(t))")
+
+function plotgif(solution::DiscSolution, tmin, tmax; kwargs...)
     ts = range(tmin, tmax, 150)
+    solver = @time solution.sol
     frames = Plots.@animate for t in ts
-        pd = dataplot(solution)
-        plot(
-            pd,
-            label = Printf.@sprintf("t = %1.2f", t),
-            ylims = (-0.1, 2.1),
-            legend = :topright,
-        )
+        pd = dataplot(solution, solver, t)
+        plot(pd ; kwargs...)
     end
     gif(frames, "temp.gif", fps = 10) |> display
 end
+
+function Base.show(io::IO, sol::DiscSolution, x)
+    println(io, string(typeof(sol).parameters[x]))
+end
\ No newline at end of file
diff --git a/src/test.jl b/src/test.jl
deleted file mode 100644
index fd312ba..0000000
--- a/src/test.jl
+++ /dev/null
@@ -1,9 +0,0 @@
-include("DiscSimulations.jl")
-using .DiscSimulationsjl
-
-initial_condition_sine(x, t, equation) = sinpi(x[1])
-source_zero(u, x, t, equations) = SVector(0, 0, 0)
-params = DiscSimulationsjl.Parameters(512, 0, 2π, initial_condition_sine, source_zero)
-
-solution = DiscSimulationsjl.main(params, DiscSimulationsjl.BurgersSimulation())
-DiscSimulationsjl.plotgif(solution, 0.0, 2.0)
\ No newline at end of file
diff --git a/src/younsi-2012.jl b/src/younsi-2012.jl
deleted file mode 100644
index e4f8b0b..0000000
--- a/src/younsi-2012.jl
+++ /dev/null
@@ -1,120 +0,0 @@
-module PolishDoughnut
-
-using DifferentialEquations
-using StaticArrays
-using ForwardDiff
-using Roots
-
-# define a derivative functor
-D(f) = x -> ForwardDiff.derivative(f, x)
-
-# metric components for the Kerr spacetime
-Σ(r, θ, a) = r^2 + a^2 * cos(θ)^2
-Δ(r, M, a) = r^2 - 2 * M * r + a^2
-function metric_components(r, θ, M, a)
-    R = 2M
-    Σ₀ = Σ(r, a, θ)
-    sinθ2 = sin(θ)^2
-
-    tt = -(1 - (R * r) / Σ₀)
-    rr = Σ₀ / Δ(r, R, a)
-    θθ = Σ₀
-    ϕϕ = sinθ2 * (r * r + a * a + (sinθ2 * R * r * a * a) / Σ₀)
-
-    tϕ = (-R * r * a * sinθ2) / Σ₀
-    @SVector [tt, rr, θθ, ϕϕ, tϕ]
-end
-
-# younsi et al. 2012 equations
-# eq. (32)
-Ω(z, M, a, n, rₖ) = (√M / (z^(3 / 2) + a * √M)) * (rₖ / z)^n
-Ω(r, θ, M, a, n, rₖ) = Ω((r * sin(θ)), M, a, n, rₖ)
-
-# eq. (30)
-Ψ₁(r, θ, M, a, Ω, Σ) = M * ((Σ - 2r^2) / (Σ^2)) * (Ω - a * sin(θ))^2 + r * sin(θ)^2
-Ψ₁(r, θ, M, a, Ω) = Ψ₁(r, θ, M, a, Ω, Σ(r, θ, a))
-
-# eq. (31)
-Ψ₂(r, θ, M, a, Ω, Σ, Δ) = sin(2θ) * ((M * r / (Σ^2)) * (a * Ω - (r^2 + a^2))^2 + Δ / 2)
-Ψ₂(r, θ, M, a, Ω) = Ψ₂(r, θ, M, a, Ω, Σ(r, θ, a), Δ(r, M, a))
-
-# eq. (28) and (29)
-function drdθ(u, p, λ)
-    # unpack parameters
-    r, θ = u
-
-    Ω₀ = inv(Ω(r, θ, p.M, p.a, p.n, p.rₖ))
-    ψ₁ = Ψ₁(r, θ, p.M, p.a, Ω₀)
-    ψ₂ = Ψ₂(r, θ, p.M, p.a, Ω₀)
-    Δ₀ = Δ(r, p.M, p.a)
-    Σ₀ = Σ(r, θ, p.a)
-
-    denom = √(Δ₀ * ψ₁^2 + ψ₂^2) * √inv(Δ₀ / Σ₀)
-    dr = ψ₂ / denom
-    dθ = -ψ₁ / denom
-
-    # use a static vector to avoid heap allocation
-    @SVector [dr, dθ]
-end
-
-function energy(r, θ, M, a, n, rₖ)
-    Ω₀ = Ω(r, θ, M, a, n, rₖ)
-    # analytic solution for energy of a circular orbit
-    g = metric_components(r, θ, M, a)
-    -(g[1] + g[5] * Ω₀) / √(-g[1] - 2g[5] * Ω₀ - g[4] * Ω₀^2)
-end
-
-function energy_minima(θ, M, a, n, rₖ, init_r)
-    # capture all but radius in a closure
-    closure = r -> energy(r, θ, M, a, n, rₖ)
-    # we want to find r for which dE/dr == 0
-    dEdr = D(closure)
-    # we give the function to find the zero for and the gradient of the function
-    # to the root finder, so we can use something like the Newton solver, which
-    # is pretty fast and accurate
-    R = find_zero((dEdr, D(dEdr)), init_r, Roots.Newton())
-    R
-end
-
-const terminate_if_negative =
-    DiscreteCallback((u, t, integrator) -> u[1] * cos(u[2]) < 0, terminate!)
-
-"""
-    isobar(M, a, n, rₖ; θ = π/2, init_r = 5.0, max_λ = 40.0)
-
-Calculate the isobar surface of a pressure supported accretion disc as in Younsi et al. (2012).
-- `M`: black hole mass,
-- `a`: (dimensionless) black hole spin,
-- `n`: angular velocity profile index of the torus,
-- `rₖ`: central Keplerian radius of the torus.
-
-The inclination of the disc from the spin axis may be set by altering `θ`. The minima of ``E(r)`` as a function of ``r``
-is numerically solved using a root finder, starting at `init_r`, and `max_λ` is the maximum integration time for the
-differential equation solver.
-"""
-function isobar(M, a, n, rₖ; θ = π / 2, init_r = 5.0, max_λ = 40.0)
-    r = energy_minima(θ, M, a, n, rₖ, init_r)
-    u0 = @SVector [r, θ]
-
-    λ_span = (0.0, max_λ)
-    # named tuple syntax
-    params = (; M = M, a = a, n = n, rₖ = rₖ)
-    prob = ODEProblem{false}(drdθ, u0, λ_span, params)
-
-    sol = solve(prob, Tsit5(), callback = terminate_if_negative, dtmax = 5e-2)
-
-    # quick and dirty unpack
-    r = first.(sol.u)
-    θ = last.(sol.u)
-
-    # then translate to cartesian coordinates
-    z = @. r * cos(θ)
-
-    # and only return those with z > 0
-    I = @. z > 0
-    r[I], z[I]
-end
-
-end # module PolishDoughnut
-
-export PolishDoughnut
diff --git a/temp.gif b/temp.gif
index 5e8fcd0..d13d63b 100644
Binary files a/temp.gif and b/temp.gif differ
diff --git a/test.jl b/test.jl
new file mode 100644
index 0000000..ae5a749
--- /dev/null
+++ b/test.jl
@@ -0,0 +1,19 @@
+include("DiscSimulations.jl")
+using .DiscSimulations
+
+initial_condition_sine(x, t, equation) = sinpi(x[1])
+#at the moment, the burgers trixi sim is expecting a function that only takes x - this could be modified for more flexibility
+initial_condition_burgers(x) = sinpi(x[1])
+
+t_span = (0.0, 2.0)
+
+source_zero(u, x, t, equations) = SVector(0, 0, 0)
+params = DiscSimulations.Parameters(512, 0.0, 2π, initial_condition_burgers, source_zero, t_span)
+
+solution = DiscSimulations.main(params, DiscSimulations.BurgersSimulation())
+s = DiscSimulations.DiscSolution(solution.sol, solution.semi, DiscSimulations.OneDimension())
+
+#NOTE: SHOW FUNCTION STILL NOT QUITE WORKING
+#DiscSimulations.Base.show(IO, solution)
+
+DiscSimulations.plotgif(s, t_span[1], t_span[2]; ylims = (-1, 1.1), legend = :topright)
\ No newline at end of file