diff --git a/FloryHugg_DiffUnbleached.ipynb b/FloryHugg_DiffUnbleached.ipynb
index 021188e63c0b784dc4a63dca3e5e923ae26e83b7..fec750d16e8b92df4185bba9d7533cfaa918af77 100644
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -66,7 +66,8 @@
     "    return c0\n",
     "\n",
     "def p_tot(p_i, p_o):\n",
-    "    return str(p_i-p_o)+'*(-0.5*tanh(3500*(x[0]-0.1))+0.5)+'+str(p_o)\n",
+    "    return str(p_i-p_o)+'*(-0.5*tanh(35000*(x[0]-0.1))+0.5)+'+str(p_o)\n",
+    "#     return '(x[0]<0.1 ? '+str(p_i)+':'+ str(p_o)+')'\n",
     "\n",
     "def create_func(f_space, expr_str, deg):\n",
     "    f = df.Function(f_space)\n",
@@ -74,7 +75,7 @@
     "    return f\n",
     "\n",
     "def eval_func(func, x):\n",
-    "    return [func([x]) for x in x]"
+    "    return np.array([func([x]) for x in x])"
    ]
   },
   {
@@ -107,6 +108,17 @@
     "plt.xlim(0, 0.2)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Set B=1\n",
+    "pi = 0.8\n",
+    "po = 0.5 - np.sqrt(0.25 + (pi**2-pi))"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -115,23 +127,34 @@
    "source": [
     "p1_i = 0.9; p1_o = 0.1\n",
     "p2_i = 0.8; p2_o = 0.2\n",
-    "p3_i = 0.9; p3_o = 0.9\n",
+    "p3_i = 0.7; p3_o = 0.3\n",
+    "p4_i = 0.9; p4_o = 0.9\n",
     "\n",
     "ct_1 = create_func(F, p_tot(p1_i, p1_o), 1)\n",
     "ct_2 = create_func(F, p_tot(p2_i, p2_o), 1)\n",
     "ct_3 = create_func(F, p_tot(p3_i, p3_o), 1)\n",
+    "ct_4 = create_func(F, p_tot(p4_i, p4_o), 1)\n",
     "\n",
     "g_1= create_func(F, '1', 1)\n",
-    "g_2= create_func(F, p_tot(p1_o/p2_o, 1), 1)\n",
-    "g_3= create_func(F, p_tot(1, 1), 1)\n",
-    "\n",
-    "c0_1 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p1_i, p1_o), 1)\n",
-    "c0_2 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p2_i, p2_o), 1)\n",
-    "c0_3 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(p3_i, p3_o), 1)\n",
+    "g_2= create_func(F, '0.5', 1)\n",
+    "g_3= create_func(F, '0.33333333333333333333', 1)\n",
+    "g_4= create_func(F, '1', 1)\n",
     "\n",
-    "c0_1 = calc_sim(c0_1, ct_1, g_1, 0.001, 50, 0)\n",
-    "c0_2 = calc_sim(c0_2, ct_2, g_2, 0.001, 50, 0)\n",
-    "c0_3 = calc_sim(c0_3, ct_3, g_3, 0.001, 50, 0)"
+    "c0_1 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p1_i, p1_o), 1)\n",
+    "c0_2 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p2_i, p2_o), 1)\n",
+    "c0_3 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p3_i, p3_o), 1)\n",
+    "c0_4 = create_func(F, 'x[0]<0.081 ? 0 :'+p_tot(p4_i, p4_o), 1)\n",
+    "for i in range(10):\n",
+    "    c0_1 = calc_sim(c0_1, ct_1, g_1, 0.001, 2, 2)\n",
+    "#     c0_2 = calc_sim(c0_2, ct_2, g_2, 0.01, 2, 2)\n",
+    "#     c0_3 = calc_sim(c0_3, ct_3, g_3, 0.01, 2, 2)\n",
+    "    c0_4 = calc_sim(c0_4, ct_4, g_4, 0.001, 2, 2)\n",
+    "    plt.plot(x, eval_func(c0_1, x), 'r')#/eval_func(c0_1, [0.099])\n",
+    "#     plt.plot(x, eval_func(c0_2, x), 'g')\n",
+    "#     plt.plot(x, eval_func(c0_3, x)/eval_func(c0_3, [0.099]), 'b')\n",
+    "    plt.plot(x, eval_func(c0_4, x), 'k')\n",
+    "    plt.xlim(0.0, 0.625)\n",
+    "    plt.ylim(0, 1.1)"
    ]
   },
   {
@@ -140,11 +163,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plt.plot(x, eval_func(c0_1, x)/np.max(eval_func(c0_1, x)))\n",
-    "plt.plot(x, eval_func(c0_2, x)/np.max(eval_func(c0_2, x)))\n",
-    "plt.plot(x, eval_func(c0_3, x)/np.array(0.47))\n",
+    "plt.plot(x, eval_func(c0_1, x)/0.9)\n",
+    "plt.plot(x, eval_func(c0_2, x)/0.8)\n",
+    "plt.plot(x, eval_func(c0_3, x)/0.7)\n",
     "plt.xlim(0.0, 0.125)\n",
-    "plt.ylim(0, 1.2)"
+    "plt.ylim(0, 1.1)"
    ]
   },
   {