diff --git a/ch01/BlackBody.ipynb b/ch01/BlackBody.ipynb
index e1b9faff..27b13e5b 100644
--- a/ch01/BlackBody.ipynb
+++ b/ch01/BlackBody.ipynb
@@ -19,12 +19,18 @@
     "from scipy.constants import h, c, k, e, N_A\n",
     "\n",
     "%matplotlib inline\n",
-    "%config InlineBackend.figure_format='retina'"
+    "%config InlineBackend.figure_format='retina'\n",
+    "\n",
+    "# Lets pick some beautiful colors\n",
+    "from matplotlib import cm\n",
+    "from cycler import cycler\n",
+    "magma = cm.get_cmap('magma', 10)  # 10 discrete colors, feel free to change\n",
+    "plt.rcParams['axes.prop_cycle'] = cycler(color=[magma(i) for i in range(magma.N)])"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 57,
+   "execution_count": null,
    "metadata": {
     "slideshow": {
      "slide_type": ""
@@ -44,24 +50,37 @@
     "    Returns:\n",
     "    - (float): Spectral radiance in W/(m^2*sr*m).\n",
     "    \"\"\"\n",
-    "    return (2.0*h*c**2) / (wavelength**5 * (np.exp((h*c)/(wavelength*k*T)) - 1))\n",
+    "    # To avoid overflow, limit the exponent to a reasonable maximum\n",
+    "    exponent = (h*c) / (wavelength*k*T)\n",
+    "    exponent = np.clip(exponent, None, 700)  # Clip at 700 to avoid overflow\n",
+    "\n",
+    "    return (2.0*h*c**2) / (wavelength**5 * (np.exp(exponent) - 1))\n",
+    "    \n",
     "\n",
-    "def plot_black_body(T=5800):\n",
+    "def plot_black_body(T, show=False):\n",
     "    \"\"\"\n",
     "    Plot the black body radiation spectrum for a given temperature T.\n",
     "    \"\"\"\n",
-    "    wavelengths = np.linspace(1e-9, 3e-6, 1000)  # Wavelength range from 1 nm to 3 um\n",
+    "    wavelengths = np.linspace(1e-9, 5e-6, 1000)  # Wavelength range from 1 nm to 3 um\n",
     "    intensities = planck(wavelengths, T)\n",
-    "\n",
-    "    plt.plot(wavelengths*1e9, intensities, label=f'T={T}')  # Convert wavelength to nm for plotting\n",
+    "    \n",
+    "    plt.plot(wavelengths*1e9, intensities, label=f'T={T}', lw=3)  # Convert wavelength to nm for plotting\n",
+    "    \n",
     "    plt.xlabel('Wavelength (nm)')\n",
-    "    plt.ylabel('Intensity (W/m^2/sr/m)')\n",
-    "    plt.grid(True)"
+    "    plt.ylabel('$Intensity (W/m^2)$')\n",
+    "    plt.grid(True)\n",
+    "    plt.legend()\n",
+    "\n",
+    "    max_int = max(planck(wavelengths, 5000)) # lets limit y axis to 5000\n",
+    "    plt.ylim(0, max_int)\n",
+    "\n",
+    "    if show: # show is needed for interactive plots\n",
+    "        plt.show() "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 53,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -99,35 +118,8 @@
     }
    ],
    "source": [
-    "plot_black_body(T=4000)\n",
-    "plot_black_body(T=5000)\n",
-    "plot_black_body(T=6000)\n",
-    "plt.legend()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 56,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "4b905a892c94445fa7d2cb4d002288fb",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "interactive(children=(IntSlider(value=5800, description='T', max=10000, min=1000, step=100), Output()), _dom_c…"
-      ]
-     },
-     "execution_count": 56,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "interactive(plot_black_body, T=(1000, 10000, 100))"
+    "for T in range(3000, 5000, 200):\n",
+    "    plot_black_body(T)"
    ]
   },
   {
@@ -135,7 +127,9 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "interactive(plot_black_body, T=(1000, 5000, 100), show=True)"
+   ]
   }
  ],
  "metadata": {