diff --git a/src/scripts/2x2_analysis.py b/src/scripts/2x2_analysis.py
new file mode 100644
index 0000000..07d26a6
--- /dev/null
+++ b/src/scripts/2x2_analysis.py
@@ -0,0 +1,48 @@
+import matplotlib.pyplot as plt
+
+
+def plot_from_file():
+    figure1, ax1 = plt.subplots()
+    figure2, ax2 = plt.subplots()
+    figure3, ax3 = plt.subplots()
+    figure4, ax4 = plt.subplots()
+
+    label = "something"
+    with open("output/test_2x2.txt") as f:
+        lines = f.readlines()
+
+        t = []
+        energy = []
+        magnetization = []
+        CV = []
+        X = []
+        for line in lines:
+            items = line.strip().split(",")
+            t.append(int(items[0]))
+            energy.append(float(items[1]))
+            magnetization.append(float(items[2]))
+            CV.append(float(items[3]))
+            X.append(float(items[4]))
+
+        ax1.plot(t, energy, label=f"<epsilon> {label}")
+        ax2.plot(t, magnetization, label=f"<|m|> {label}")
+        ax3.plot(t, CV, label=f"CV {label}")
+        ax4.plot(t, X, label=f"X {label}")
+
+        # ax1.set_yscale("log")
+        # ax2.set_yscale("log")
+
+    figure1.legend()
+    figure1.savefig("../latex/images/2x2_energy.pdf")
+    figure2.legend()
+    figure2.savefig("../latex/images/2x2_magnetization.pdf")
+    figure3.legend()
+    figure3.savefig("../latex/images/2x2_heat_capacity.pdf")
+    figure4.legend()
+    figure4.savefig("../latex/images/2x2_susceptibility.pdf")
+
+def main():
+    plot_from_file()
+
+if __name__ == "__main__":
+    main()
diff --git a/src/scripts/burn_in_time.py b/src/scripts/burn_in_time.py
new file mode 100644
index 0000000..cd86a97
--- /dev/null
+++ b/src/scripts/burn_in_time.py
@@ -0,0 +1,40 @@
+import matplotlib.pyplot as plt
+
+def plot_from_file():
+    files = [
+        "output/burn_in_time_1_0.txt",
+        "output/burn_in_time_2_4.txt",
+    ]
+    labels = [
+        "1.0",
+        "2.4"
+    ]
+    figure1, ax1 = plt.subplots()
+    figure2, ax2 = plt.subplots()
+
+    for infile, label in zip(files, labels):
+        with open(infile) as f:
+            lines = f.readlines()
+
+            t = []
+            energy = []
+            magnetization = []
+            for line in lines:
+                items = line.strip().split(",")
+                t.append(int(items[0]))
+                energy.append(float(items[1]))
+                magnetization.append(float(items[5]))
+
+            ax1.plot(t, energy, label=f"<epsilon> {label}")
+            ax2.plot(t, magnetization, label=f"<|m|> {label}")
+
+    figure1.legend()
+    figure1.savefig("../latex/images/burn_in_time_energy.pdf")
+    figure2.legend()
+    figure2.savefig("../latex/images/burn_in_time_magnetization.pdf")
+
+def main():
+    plot_from_file()
+
+if __name__ == "__main__":
+    main()
diff --git a/src/scripts/pd_estimate.py b/src/scripts/pd_estimate.py
new file mode 100644
index 0000000..c5cd327
--- /dev/null
+++ b/src/scripts/pd_estimate.py
@@ -0,0 +1,24 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from pathlib import Path
+
+def plot(infile, outfile):
+
+    figure1, ax1 = plt.subplots()
+    arr = []
+    with open(infile) as f:
+        lines = f.readlines()
+        for line in lines:
+            vals = line.strip().split(",")
+            arr.append(float(vals[0]))
+
+    ax1.hist(arr, np.arange(min(arr), max(arr) + .02, .02), density=True)
+
+    figure1.savefig(outfile)
+
+if __name__ == "__main__":
+    plot("output/pd_estimate/estimate_1_0.txt", "../latex/images/pd_estimate_1_0.pdf")
+    plot("output/pd_estimate/estimate_2_4.txt", "../latex/images/pd_estimate_2_4.pdf")
+
+
+
diff --git a/src/scripts/phase_transition.py b/src/scripts/phase_transition.py
new file mode 100644
index 0000000..7782041
--- /dev/null
+++ b/src/scripts/phase_transition.py
@@ -0,0 +1,61 @@
+import matplotlib.pyplot as plt
+
+from pathlib import Path
+
+def plot_phase_transition(indir, outdir):
+    files = [
+        "size_20.txt",
+        "size_40.txt",
+        "size_60.txt",
+        "size_80.txt",
+        "size_100.txt",
+    ]
+    labels = [
+        "L = 20",
+        "L = 40",
+        "L = 60",
+        "L = 80",
+        "L = 100",
+    ]
+
+    figure1, ax1 = plt.subplots()
+    figure2, ax2 = plt.subplots()
+    figure3, ax3 = plt.subplots()
+    figure4, ax4 = plt.subplots()
+
+    for file, label in zip(files, labels):
+        t = []
+        e = []
+        m = []
+        CV = []
+        X = []
+
+        with open(Path(indir, file)) as f:
+            lines = f.readlines()
+            for line in lines:
+                l = line.strip().split(",")
+                t.append(float(l[0]))
+                e.append(float(l[1]))
+                m.append(float(l[2]))
+                CV.append(float(l[3]))
+                X.append(float(l[4]))
+
+        ax1.plot(t,e,label=label)
+        ax2.plot(t,m,label=label)
+        ax3.plot(t,CV, label=label)
+        ax4.plot(t,X, label=label)
+
+    figure1.legend()
+    figure2.legend()
+    figure3.legend()
+    figure4.legend()
+
+    figure1.savefig(Path(outdir, "energy.pdf"))
+    figure2.savefig(Path(outdir, "magnetization.pdf"))
+    figure3.savefig(Path(outdir, "heat_capacity.pdf"))
+    figure4.savefig(Path(outdir, "susceptibility.pdf"))
+
+if __name__ == "__main__":
+    plot_phase_transition("fox_output/phase_transition/10M/", "../latex/images/phase_transition/wide/10M/")
+    plot_phase_transition("fox_output/phase_transition/1M/", "../latex/images/phase_transition/wide/1M/")
+    plot_phase_transition("fox_output/phase_transition/narrow/10M/", "../latex/images/phase_transition/narrow/10M/")