diff --git a/Final/Final - Tight-binding propagation method.ipynb b/Final/Final - Tight-binding propagation method.ipynb
index dc3af21..520d5e0 100644
--- a/Final/Final - Tight-binding propagation method.ipynb	
+++ b/Final/Final - Tight-binding propagation method.ipynb	
@@ -308,8 +308,8 @@
    "source": [
     "def integrate(yk, x):\n",
     "    \"\"\"\n",
-    "    Numerically integrates function yk over [x[0], x[-1]] using\n",
-    "    Simpson's rule over the grid provided by x.\n",
+    "    Numerically integrates function yk over [x[0], x[-1]] using Simpson's 3/8\n",
+    "    rule over the grid provided by x.\n",
     "    \n",
     "    Args:\n",
     "        yk: function of one numerical argument that returns a numeric\n",
@@ -326,23 +326,17 @@
     "    if callable(yk):\n",
     "        yk = yk(x)\n",
     "    \n",
-    "    #a, b = yk[0], yk[-1]\n",
-    "    \n",
+    "    # The distance h_i = x[i + 1] - x[i] is not necessarily constant. The choice of\n",
+    "    # partitioning of the interval is subject to mathematical considerations I will\n",
+    "    # not go into.\n",
     "    h = x[1:] - x[:-1]\n",
     "    \n",
-    "    # Instead of summing over something with n/2, we use step size 2,\n",
-    "    # thus avoiding any issues with 2 ∤ n.\n",
-    "    #integral = h/3*(a + 2*np.sum(yk[2:-1:2]) + 4*np.sum(yk[1:-1:2]) + b)\n",
-    "    \n",
-    "    #integral = 3*h/8*(a + 3*np.sum(yk[2:-1:2]) + 3*np.sum(yk[1:-1:2]) + b)\n",
-    "    #integral = (b - a)/8*(x[0] + x[-1]) + 3*(x[1:] - x[:-1])/8*(  )\n",
     "    integral = 0\n",
     "    integral += 3/8*(x[1] - x[0])*yk[0]\n",
-    "    integral += 9/8*yk[1:-1:3]@h[1::3]\n",
-    "    integral += 9/8*yk[2:-1:3]@h[2::3]\n",
-    "    integral += 6/8*yk[:-1:3]@h[::3]\n",
+    "    integral += 9/8*h[1::3]@yk[1:-1:3]\n",
+    "    integral += 9/8*h[2::3]@yk[2:-1:3]\n",
+    "    integral += 6/8*h[ ::3]@yk[ :-1:3]\n",
     "    integral += 3/8*(x[-1] - x[-2])*yk[-1]\n",
-    "    #integral = 3/8*( (x[1] - x[0])*yk[0] + np.sum(yk[]) + (x[-1] - x[-2])*yk[-1])\n",
     "    return integral"
    ]
   },
@@ -363,22 +357,20 @@
      "task": false
     }
    },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Exact: 333333333.3333333\n",
-      "Integrated: 333333448.1291072\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "def test_integrate():\n",
+    "    # Test integral 1 of f with F its primitive with integration constant 0\n",
     "    f = lambda x: x**2\n",
-    "    x = np.logspace(0, 3, 10000000)\n",
-    "    print(\"Exact:\", 1000**3/3)\n",
-    "    print(\"Integrated:\", integrate(f, x))\n",
+    "    F = lambda x: x**3/3\n",
+    "    x = np.logspace(0, 3, 1000000)\n",
+    "    assert np.isclose(integrate(f, x), F(x[-1]) - F(x[0]))\n",
+    "    \n",
+    "    # Test integral 2 of f with F its primitive with integration constant 0\n",
+    "    f = lambda x: np.sin(2*x)/(2 + np.cos(2*x))\n",
+    "    F = lambda x: -.5*np.log(np.cos(2*x) + 2)\n",
+    "    x = np.linspace(0, 10, 1000)\n",
+    "    assert np.isclose(integrate(f, x), F(x[-1]) - F(x[0]))\n",
     "    \n",
     "test_integrate()"
    ]