diff --git a/Fig2a_plot.py b/Fig2a_plot.py
new file mode 100644
index 0000000..bd66e37
--- /dev/null
+++ b/Fig2a_plot.py
@@ -0,0 +1,389 @@
+def onorm(M):
+    MM=np.conj(M).T@M
+    eigenvalues,eigenvectors=np.linalg.eig(MM)
+    
+    return np.sqrt(np.max(eigenvalues))
+
+
+
+
+
+
+
+import math as math
+import numpy as np
+import sympy as sp
+from mpmath import coth
+import scipy.linalg as scln
+from scipy.integrate import trapz
+import itertools
+from sympy.tensor.array.expressions import ArraySymbol
+from sympy.abc import i, j, k
+import h5py
+import matplotlib.pyplot as plt
+from engine.inverse.inverse import *
+plt.rc("font", family="serif")
+plt.rc("xtick")
+
+
+plt.rc("ytick")
+fig = plt.figure(  )  # width=8 inches, height=6 inches
+
+ax = fig.add_subplot(1, 1, 1)
+N=500
+wantdelt=0.1
+t=np.linspace(0,wantdelt*(N-1),num=N)
+
+Idval=np.array([[1,0],[0,1]],dtype = 'complex')
+Hs=np.array([[0,1],[1,0]],dtype = 'complex')
+delt=wantdelt
+G=scln.expm(np.kron(1j*Hs*delt/2,Idval)+np.kron(Idval,(-1j*Hs*delt/2)))
+filename = f'Data/Fig 2 Data/jobgqme0.11.0/dyn0.1beta5delt0.1ohmicity1.0.h5'
+
+h5 = h5py.File(filename,'r')
+kmax=14
+#MK=np.zeros((4,4,kmax+1),dtype='complex')
+M = h5["M"]  # VSTOXX futures data
+A = np.array(h5["10AT"])
+print(list(A))
+UU=np.zeros((4,4,N),dtype = 'complex')
+UU[:,:,1]=M[:,:,1]
+Mnorm=np.zeros((kmax+1),dtype = 'complex')
+MMM=np.zeros((4,4,kmax+1),dtype = 'complex')
+
+for k in range(2,kmax+1):
+    tempM=np.zeros((4,4),dtype = 'complex')
+    for j in range(1,k):
+        tempM += np.matmul(M[:, :, j],UU[:, :, k - j])
+    (UU[:,:,(k)])=M[:,:,(k)]+tempM
+for i in range(1,kmax+1):
+    UU[:,:,i]=np.matmul(G,np.matmul(UU[:,:,i],G))
+    MMM[:,:,i]=M[:,:,i]
+MK=KfromU(UU,kmax)
+#ax.plot(np.arange(0,14)*wantdelt,MK[1,2,:],color="k", ls="solid",label="Dyck/INFPI $\mathcal{K}_{12}^N$")
+#ax.plot(np.arange(0,14)*wantdelt,MK[0,0,:],color=(0.6,0.6,0.6), ls="solid",label="Dyck/INFPI $\mathcal{K}_{00}^N$")
+#ax.plot(np.arange(0,14)*wantdelt,MK[0,3,:],color=(0.2,0.2,0.2), ls="solid",label="Dyck/INFPI $\mathcal{K}_{03}^N$")
+for i in range(kmax+1):
+    Mnorm[i]=onorm(M[:,:,i]/(delt**2))
+#ax.plot(t,A[0,:],color="k")
+ax.semilogy(np.arange(2,kmax+1)*wantdelt,Mnorm[2:],color="k", ls="none",marker='.',markersize=7)
+
+h5.close()
+
+ax.set_facecolor((0.9,1,1))
+plt.xlabel("Time",fontsize=15)
+#ax.set_ylabel(r"$\sigma_z$")
+#ax.set_ylim([-0.2,1])
+ax.set_ylim([0.0001,15])
+ax.set_xlim([0,4])
+
+#plt.legend()
+#plt.savefig("plots/K051.png",dpi=1000)
+#plt.show()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+import math as math
+import numpy as np
+import sympy as sp
+from mpmath import coth
+import scipy.linalg as scln
+from scipy.integrate import trapz
+import itertools
+from sympy.tensor.array.expressions import ArraySymbol
+from sympy.abc import i, j, k
+import h5py
+import matplotlib.pyplot as plt
+from engine.inverse.inverse import *
+filename = f'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.1s1initial1.h5'
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+A1 = np.transpose(h5["rho"])
+
+filename = f'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.1s1initial2.h5'
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+A2 = np.transpose(h5["rho"])
+
+filename = f'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.1s1initial3.h5'
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+A3 = np.transpose(h5["rho"])
+
+filename = 'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.1s1initial4.h5'
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+A4 = np.transpose(h5["rho"])
+
+Uex=UfromA(A1,A2,A3,A4)
+A5=A1.copy()
+for i in range(1,len(A1)):  
+    A5[:,i]=np.matmul(Uex[:,:,i],[1,0,0,0])
+mmax=150
+Mex=KfromU(Uex,mmax)
+
+t=np.arange(0,0.1*(mmax+1),0.1)
+ddelt=0.1
+Mnormex=np.zeros((mmax+1),dtype = 'complex')
+for i in range(mmax+1):
+    Mnormex[i]=onorm(Mex[:,:,i]/(ddelt**2))
+ax.plot(t[2:],Mnormex[2:],color=(0,0,0),alpha=1, ls="solid",label="HEOM $\mathcal{K}_{12}^N$")
+plt.xticks(fontsize=15)
+
+plt.yticks(fontsize=15)
+#plt.savefig("plots/K051HEOM.png",dpi=1000)
+#plt.plot(t[0:500],A1[0,:500])
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+import math as math
+import numpy as np
+import sympy as sp
+import scipy as sc
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import scipy.linalg as scln
+from scipy.integrate import trapz, quad
+import itertools
+from sympy.tensor.array.expressions import ArraySymbol
+from sympy.abc import i, j, k
+mpl.rcParams.update(mpl.rcParamsDefault)
+beta=5
+kmax=40
+N=500
+wantdelt=0.1
+kondo=0.1
+ga=kondo* np.pi / 2
+No=20000
+hbar=1
+om=np.linspace(-100,100,No)
+yy1=np.zeros(No,dtype = "complex")
+yy2=np.zeros(No,dtype = "complex")
+tarr = np.linspace(0,wantdelt*(N-1),num=N)
+delt=(tarr[1]-tarr[0])
+Hs=np.array([[0,1],[1,0]],dtype = "complex")
+wcut = 7.5
+JJ=om.copy()
+ohmicity=1
+wc=wcut
+def J(xs):  # Spectral density (Ohmic bath with exponential cutoff)
+    output = np.zeros_like(xs)
+    output[xs > 0] = ga * (xs**ohmicity)[xs > 0] * np.exp(-xs[xs > 0] / wc)/(wc**(ohmicity-1))
+    output[xs <= 0] = - ga * ((-xs)**ohmicity)[xs <= 0] * np.exp(xs[xs <= 0] / wc)/(wc**(ohmicity-1))
+    #output=np.heaviside(xs,0.5)ga*(xs**ohmicity)/(wc**(1-ohmicity))*np.exp(-(xs)/wc)
+    return output
+for i in enumerate(om):
+    JJ[i[0]]=J(i[1])
+
+etakkd=np.zeros(kmax+1,dtype="complex")
+
+
+delt2=delt*1
+for iii in enumerate(om):   
+    yy1[iii[0]]=1/(2*np.pi)*J(om[iii[0]])/(om[iii[0]]**2)*np.exp(beta*hbar*om[iii[0]]/2)/(math.sinh(beta*hbar*om[iii[0]]/2))*(1-np.exp(-1j*om[iii[0]]*delt2))
+etanul=trapz(yy1,om)
+
+
+onetwo=[1,-1]
+def yy2f(x,diffk): 
+    return 2/(1*np.pi)*J(x)/(x**2)*np.exp(beta*hbar*x/2)/(math.sinh(beta*hbar*x/2))*((math.sin(x*delt2/2))**2)*np.exp(-1j*x*delt2*(diffk))
+for difk in np.linspace(1,kmax,kmax,dtype="int"):
+    for i in enumerate(om):
+        #print(om[i[0]])
+        yy2[i[0]]=2/(1*np.pi)*J(om[i[0]])/(om[i[0]]**2)*np.exp(beta*hbar*om[i[0]]/2)/(math.sinh(beta*hbar*om[i[0]]/2))*((math.sin(om[i[0]]*delt2/2))**2)*np.exp(-1j*om[i[0]]*delt2*(difk))
+    #aaa=quad(yy2f,-np.inf,np.inf,args=difk)
+    etakkd[difk]=trapz(yy2,om)
+def Influencediff(dx1,dx2,dx3,dx4,diff):
+    x1=onetwo[dx1]
+    x2=onetwo[dx2]
+    x3=onetwo[dx3]
+    x4=onetwo[dx4]
+    Sum=-1/hbar*(x3-x4)*(etakkd[diff]*x1-np.conjugate(etakkd[diff])*x2)     # eq 12 line 1 Nancy quapi I           
+    return np.exp(Sum)
+
+def Influencenull(dx1,dx2): # eq 12 line 2 Nancy quapi I    
+    x1=onetwo[dx1]
+    
+    x2=onetwo[dx2]
+    Sum=-1/hbar*(x1-x2)*((etanul)*x1-np.conjugate(etanul)*x2)                
+    return np.exp(Sum)
+    
+
+I0_val=np.array([Influencenull(0,0),Influencenull(0,1),Influencenull(1,0),Influencenull(1,1)], dtype = "complex")
+
+TI_val=np.zeros((4,4,kmax+1),dtype = "complex")  
+I_val=np.zeros((4,4,kmax+1),dtype = "complex")  
+A=np.zeros((4,len(tarr)),dtype = "complex")
+A[0,0]=1
+P_val=np.zeros((4,4),dtype = "complex")
+
+
+for i in np.arange(0,kmax+1):
+    I_val[0,0,i]=Influencediff(0,0,0,0,i)
+    I_val[0,1,i]=Influencediff(1,0,0,0,i)
+    I_val[0,2,i]=Influencediff(0,1,0,0,i)
+    I_val[0,3,i]=Influencediff(1,1,0,0,i)
+
+    I_val[1,0,i]=Influencediff(0,0,1,0,i)
+    I_val[1,1,i]=Influencediff(1,0,1,0,i)
+    I_val[1,2,i]=Influencediff(0,1,1,0,i)
+    I_val[1,3,i]=Influencediff(1,1,1,0,i)
+
+    I_val[2,0,i]=Influencediff(0,0,0,1,i)
+    I_val[2,1,i]=Influencediff(1,0,0,1,i)
+    I_val[2,2,i]=Influencediff(0,1,0,1,i)
+    I_val[2,3,i]=Influencediff(1,1,0,1,i)
+
+    I_val[3,0,i]=Influencediff(0,0,1,1,i)
+    I_val[3,1,i]=Influencediff(1,0,1,1,i)
+    I_val[3,2,i]=Influencediff(0,1,1,1,i)
+    I_val[3,3,i]=Influencediff(1,1,1,1,i)
+for i in range(kmax+1):
+    TI_val[:,:,i]=I_val[:,:,i]-np.ones(4)
+#plt.rc('text', usetex=True)
+
+#plt.rc("text", usetex=True)
+
+
+
+I0_valm=np.zeros((4,4),dtype="complex")
+for v in range(4):
+    I0_valm[v,v]=I0_val[v]  
+#fnorm[0]=np.linalg.norm(I_val)
+print(np.linalg.norm(I0_valm))
+x=np.log(np.arange(1,kmax+1))
+for i in range(kmax):
+    x[i]=onorm(TI_val[:,:,i])
+#print(sc.stats.linregress(x, y=np.log(fnorm[1:]), alternative="two-sided"))
+#ax.semilogy(range(1,kmax+1)*delt,np.abs(I_val[1,2,1:]-1),color="k", ls="solid",label="Re${I}_{12}^{\Delta k}$")
+#ax.semilogy(range(1,kmax+1)*delt,np.abs(I_val[1,3,1:]-1),color=(0.6,0.6,0.6), ls="solid",label="Re${I}_{13}^{\Delta k}$")
+#ax.semilogy(range(1,kmax+1)*delt,np.abs(I_val[0,0,1:]-1),color=(0.2,0.2,0.2), ls="solid",label="Re${I}_{00}^{\Delta k}$")
+
+#ax.semilogy(range(1,14)*delt,np.abs(I_val[1,2,1:14]-1),color="k", ls="none",marker='x')
+#ax.semilogy(range(1,14)*delt,np.abs(I_val[1,3,1:14]-1),color=(0.6,0.6,0.6), ls="none",marker='x')
+#ax.semilogy(range(1,14)*delt,np.abs(I_val[0,0,1:14]-1),color=(0.2,0.2,0.2), ls="none",marker='x')
+
+#ax.set_xlabel("Time")
+#ax.set_xlim([delt,delt*kmax])
+#ax.set_ylim([0.99,1.01])
+#plt.plot(np.log(abs(etakkd[1:])))
+#plt.plot(2.5*1/x**0.5
+#ax.set_ylabel(r'$\|\mathcal{K}_{N \Delta t}\|_O$', color='k',fontsize=15)
+#ax2 = ax.twinx()
+ax.semilogy(range(2,kmax+1)*delt,x[1:], color='#5E091A')
+ax.semilogy(range(2,14+1)*delt,x[1:14], color='#5E091A', ls="none",marker='x',markersize=7)
+#ax.set_ylabel(r'$\|\tilde{I}\|_O$', color='#5E091A',fontsize=15)
+ax.tick_params('y', color='#5E091A')
+#fnorm=np.zeros((kmax+1),dtype="complex")
+#matrixone=np.ones((4),dtype="complex")
+#plt.legend()
+
+plt.xticks(fontsize=15)
+plt.yticks(fontsize=15)
+filename="Data/Fig 2 Data/crestcrest6/dyn0.15707963267948966ohmicity1.0.h5"
+
+h5 = h5py.File(filename,'r')
+#kmax=14
+sett = h5["sett"]  # VSTOXX futures data
+nnn = h5["norm"]
+nnnn=np.zeros(17,dtype='complex')
+aaa=np.zeros(17,dtype='complex')
+bbb=np.zeros((4,4,18),dtype='complex')
+#plt.plot(range(2,18)*delt,list(sett[1:]),alpha=0.9,color='black', ls="dashed",label=r'$\|\mathcal{K}_{N,crest}\|$')
+for i in range(17):
+    
+    bbb[:,:,i]=sett[:,:,i]
+    #aaa[i]=onorm(bbb[:,:,i]-MMM[:,:,i+1]/0.01)
+    nnnn[i]=nnn[i]
+h5.close()
+print(Mex[2,1,1:8]/0.01)
+print(bbb[2,1,:7])
+
+plt.plot(range(2,18)*delt,(nnnn[1:]),alpha=0.9,color='black', ls="dashed",label=r'$\|\mathcal{K}_{N,crest}-\mathcal{K}_{N}\|$')
+#plt.plot(range(2,16)*delt,(aaa[0:-3]),alpha=0.9,color='black', ls="dashed",label=r'$\|\mathcal{K}_{N,crest}-\mathcal{K}_{N}\|$')
+plt.savefig(f"Fig1a.eps",format='eps',bbox_inches="tight")
+
+plt.show()
\ No newline at end of file
diff --git a/Fig2b_plot.py b/Fig2b_plot.py
index cb6629f..8d02bbc 100644
--- a/Fig2b_plot.py
+++ b/Fig2b_plot.py
@@ -17,7 +17,7 @@ def onorm(M):
 from sympy.abc import i, j, k
 import h5py
 import matplotlib.pyplot as plt
-from pyeom.inverse.inverse import *
+from engine.inverse.inverse import *
 plt.rc("font", family="serif")
 plt.rc("xtick")
 
@@ -145,7 +145,7 @@ def onorm(M):
 from sympy.abc import i, j, k
 import h5py
 import matplotlib.pyplot as plt
-from pyeom.inverse.inverse import *
+from engine.inverse.inverse import *
 filename = f'Data/Fig 2 Data/heompaper10.115.0/heompaper10.115.0/tmax15.0beta5alpha0.5s1initial1.h5'
 
 h5 = h5py.File(filename,'r')
diff --git a/Fig2c_plot.py b/Fig2c_plot.py
index c58b794..84e040d 100644
--- a/Fig2c_plot.py
+++ b/Fig2c_plot.py
@@ -21,7 +21,7 @@ def onorm(M):
 from sympy.abc import i, j, k
 import h5py
 import matplotlib.pyplot as plt
-from pyeom.inverse.inverse import *
+from engine.inverse.inverse import *
 plt.rc("font", family="serif")
 plt.rc("xtick")
 
@@ -149,7 +149,7 @@ def onorm(M):
 from sympy.abc import i, j, k
 import h5py
 import matplotlib.pyplot as plt
-from pyeom.inverse.inverse import *
+from engine.inverse.inverse import *
 filename = f'Data/Fig 2 Data/ttt/ttt/heomdynamicssubohmicpapertrueyes0.110.0/tmax10.0beta5alpha0.5s0.5initial1.h5'
 
 h5 = h5py.File(filename,'r')
diff --git a/Fig3_panels_a1b1c1.py b/Fig3_panels_a1b1c1.py
index df8f54b..7b70d0f 100644
--- a/Fig3_panels_a1b1c1.py
+++ b/Fig3_panels_a1b1c1.py
@@ -11,7 +11,7 @@
 from sympy.abc import i, j, k
 import h5py
 import matplotlib.pyplot as plt
-from pyeom.inverse.inverse import *
+from engine.inverse.inverse import *
 plt.rc("font", family="serif")
 plt.rc("xtick")
 plt.rc("ytick")
@@ -169,7 +169,7 @@
 from sympy.abc import i, j, k
 import h5py
 import matplotlib.pyplot as plt
-from pyeom.inverse.inverse import *
+from engine.inverse.inverse import *
 filename = f'Data\heomforK0.015.0/tmax5.0beta5alpha0.5s0.5initial1.h5'
 
 h5 = h5py.File(filename,'r')
diff --git a/Fig3_panels_a2b2c2_theoretical.py b/Fig3_panels_a2b2c2_theoretical.py
index e5ab022..29ae08b 100644
--- a/Fig3_panels_a2b2c2_theoretical.py
+++ b/Fig3_panels_a2b2c2_theoretical.py
@@ -15,6 +15,8 @@
 from scipy.interpolate import interp1d
 ##This code is skipping the rho->U->K because they are exact in principle
 ###Set kondo and ohmicity below to plot the respective plots
+####Modify this code (and the from data code) accordingly to produce Figs. S3 and S6
+
 
 kkmax=15
 N=500
diff --git a/Fig_S4a_plot.py b/Fig_S4a_plot.py
new file mode 100644
index 0000000..c92b4bd
--- /dev/null
+++ b/Fig_S4a_plot.py
@@ -0,0 +1,100 @@
+import matplotlib.pyplot as plt
+import h5py
+import numpy as np
+file_path = 'Data\Fig S4 Data\Driven_HEOM_results\driven_heom_d_eps_1.out'
+
+data = []  # to store the parsed data
+
+with open(file_path, 'r') as file:
+    for line in file:
+        # Split each line into a list of values using space as the delimiter
+        values = [float(val) for val in line.split()]
+        data.append(values)
+
+# Now 'data' is a list of lists, where each inner list represents a line of values
+# You can access specific values using indexing, e.g., data[0] for the first line
+# or data[0][1] for the second value in the first line.
+
+# Example: Print the entire data
+time=[]
+y=[]
+y1=[]
+y2=[]
+for line in data:
+    time.append((line[0]))
+    y.append((line[1]))
+    y1.append((line[7]))
+    y2.append((line[3]))
+fig, ax = plt.subplots()
+ax.set_facecolor((0.95,0.95,0.95))
+plt.rc("font", family="serif")
+plt.plot(time,y,color='red',ls='dashed',marker='o',markevery=15)
+plt.plot(time,y1,color='red',ls='dashed',marker='x',markevery=15)
+plt.plot(time,y2,color='red',ls='dashed',marker='+',markevery=15)
+filename='Data\Fig S4 Data\drivendata/d.h5'
+h5 = h5py.File(filename,'r')
+#kmax=14
+kmax=6
+N=200
+wantdelt=0.1
+tarr = np.linspace(0,wantdelt*(N-1),num=N)
+M = h5["A"]  # VSTOXX futures data
+
+plt.plot(tarr,M[0,:],color='#808080',marker='o',ls='dotted',markevery=20,label='Dyck, order 4')
+
+plt.plot(tarr,M[3,:],color='#808080',marker='x',ls='dotted',markevery=20)
+plt.plot(tarr,M[2,:],color='#808080',marker='+',ls='dotted',markevery=20)
+
+
+
+
+
+filename='Data\Fig S4 Data\drivendata/td25.h5'
+h5 = h5py.File(filename,'r')
+#kmax=14
+kmax=4
+N=200
+wantdelt=0.1
+tarr = np.linspace(0,wantdelt*(N-1),num=N)
+M = h5["A"]  # VSTOXX futures data
+UG= h5["UB"]
+MM=np.zeros((4,200),dtype='complex')
+for i in np.arange(0,N):
+    MM[:,i]=M[:,i]
+for i in np.arange(kmax+1,N):
+    MM[:,i]=np.matmul(UG[:,:,i],np.transpose([1,0,0,0]))
+plt.plot(tarr,MM[0,:],color='#808080',marker='o',markevery=20,label='Dyck, order 4 + crest term upto 5')
+
+plt.plot(tarr,MM[3,:],color='#808080',marker='x',markevery=20)
+plt.plot(tarr,MM[2,:],color='#808080',marker='+',markevery=20)
+
+
+
+
+
+filename='Data\Fig S4 Data\drivendata/td2.h5'
+h5 = h5py.File(filename,'r')
+#kmax=14
+kmax=4
+N=200
+wantdelt=0.1
+tarr = np.linspace(0,wantdelt*(N-1),num=N)
+M = h5["A"]  # VSTOXX futures data
+UG= h5["UB"]
+MM=np.zeros((4,200),dtype='complex')
+for i in np.arange(0,N):
+    MM[:,i]=M[:,i]
+for i in np.arange(kmax+1,N):
+    MM[:,i]=np.matmul(UG[:,:,i],np.transpose([1,0,0,0]))
+plt.plot(tarr,MM[0,:],color='black',marker='o',markevery=20,label='Dyck, order 4 + crest term upto 6')
+
+plt.plot(tarr,MM[3,:],color='black',marker='x',markevery=20)
+plt.plot(tarr,MM[2,:],color='black',marker='+',markevery=20)
+plt.xticks(fontsize=15, fontname = 'serif')
+plt.yticks(fontsize=15, fontname = 'serif')
+plt.xlabel('Time',fontsize=15, fontname = 'serif')
+plt.ylim([-0.5,1])
+plt.xlim([0,20])
+plt.legend()
+#plt.savefig("driven5.png",dpi=1000)
+plt.show()
\ No newline at end of file
diff --git a/README.md b/README.md
index 7224c34..ecc7fa5 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,5 @@
 # pyeom
-A python code package accompanying "Unified Framework for Open Quantum Dynamics with Memory" at Nature Communications by Ivander, Lindoy, and Lee.
+Python-based real-time quantum dynamics methods
 
-The main directory contains the code to plot the figures in the main text (and a figure set in the supplementary information). They are named by the figures they represent in the manuscript. Some will need modifying parameters to specify which figure panel the codes will produce. Important information are commented in the early parts of the code.
 
-Inside "Data/" you will find the raw data for each figures and the job files that were required to produce them.
-
-Inside "pyeom/" you will find the engine behind the computations, with codes for performing HEOM, QuaPI, Dyck path, TTM simulations, etc.
+call Qua.py [kmax]