diff --git a/ps11-2.13-2.py b/ps11-2.13-2.py
deleted file mode 100644
index c80c494..0000000
--- a/ps11-2.13-2.py
+++ /dev/null
@@ -1,20 +0,0 @@
-import numpy as np
-import pylab as pl
-
-x_max = 10.0
-x_vals = np.linspace(-x_max, x_max, 1000)
-dx = x_vals[1] - x_vals[0]
-gaussian = np.exp(-x_vals * x_vals / 2)
-
-psi_funcs = [gaussian,
-             x_vals * gaussian,
-             (1 - 2 * x_vals * x_vals) * gaussian,
-             x_vals * (1 - (2.0 / 3.0) * x_vals * x_vals) * gaussian]
-
-for psi in psi_funcs:
-    a = 1.0 / (np.sqrt(psi.dot(psi)) * dx)
-    psi *= a
-
-for psi in psi_funcs:
-    pl.plot(x_vals, psi)
-    pl.show()
diff --git a/ps12-2.py b/ps12-2.py
new file mode 100644
index 0000000..e4991e5
--- /dev/null
+++ b/ps12-2.py
@@ -0,0 +1,56 @@
+import numpy as np
+import pylab as pl
+
+
+# recurrence relation is H_n(x) = 2*x*H_{n-1}(x) - 2*(n-1)*H_{n-2}(x)
+# source https://en.wikipedia.org/wiki/Hermite_polynomials
+def hermite(x: np.ndarray, n: int) -> np.ndarray:
+    if n == 0:
+        return 1.0 + 0.0 * x
+    if n == 1:
+        return 2.0 * x
+    return 2.0 * x * hermite(x, n - 1) - 2.0 * (n - 1) * hermite(x, n - 2)
+
+
+def normalize(psi: np.ndarray) -> None:
+    psi /= np.sqrt(np.sum(psi * np.conjugate(psi)))
+
+
+x_vals = np.linspace(-5, 5, 1000)
+gaussian = np.exp(-x_vals ** 2)
+
+psi_funcs = []
+for n in range(4):
+    psi = hermite(x_vals, n) * gaussian
+    normalize(psi)
+    psi_funcs.append(psi)
+
+print("===================")
+print("Check Orthogonality")
+print("===================")
+for i in range(4):
+    for j in range(i, 4):
+        print(f'psi_{i} * psi_{j} = {np.dot(psi_funcs[i], psi_funcs[j])}')
+
+print("==================")
+print("Expectation Values")
+print("==================")
+for psi in psi_funcs:
+    pass  # TODO
+
+pl.rcParams['figure.dpi'] = 300
+fig, axs = pl.subplots(2, 2, tight_layout=True)
+fig.tight_layout(pad=2.0)
+axs[0, 0].plot(x_vals, psi_funcs[0])
+axs[0, 0].set_title('ψ0')
+axs[0, 0].grid()
+axs[0, 1].plot(x_vals, psi_funcs[1], 'tab:orange')
+axs[0, 1].set_title('ψ1')
+axs[0, 1].grid()
+axs[1, 0].plot(x_vals, psi_funcs[2], 'tab:green')
+axs[1, 0].set_title('ψ2')
+axs[1, 0].grid()
+axs[1, 1].plot(x_vals, psi_funcs[3], 'tab:red')
+axs[1, 1].set_title('ψ3')
+axs[1, 1].grid()
+pl.show()