diff --git a/FloryHugg_DiffUnbleached.ipynb b/FloryHugg_DiffUnbleached.ipynb
index 6c77811d0bf5e36a84a9a1ab086696d681c09340..c87251276bca4db12ab6ed1f416f557d9c8822bb 100644
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -93,37 +93,40 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "bp = 0.1 # Boundary position at IC\n",
+    "sym = 2 # Symmetry of the problem\n",
+    "dt = 0.0001\n",
+    "\n",
     "c_tot_1 = create_func(F, p_tot(0.9, 0.9), 1)\n",
-    "c0_1 = create_func(F, 'x[0]<0.1 ? 0 :' + p_tot(0.9, 0.9), 1)\n",
+    "c0_1 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' + p_tot(0.9, 0.9), 1)\n",
     "Ga0_1 = create_func(F, p_tot(1, 1/9), 1)\n",
     "\n",
     "c_tot_2 = create_func(F, p_tot(0.9, 0.45), 1)\n",
-    "c0_2 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.45), 1)\n",
+    "c0_2 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.45), 1)\n",
     "Ga0_2 = create_func(F,p_tot(1, 0.08), 1)\n",
     "\n",
     "c_tot_3 = create_func(F, p_tot(0.9, 0.3), 1)\n",
-    "c0_3 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.3), 1)\n",
+    "c0_3 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.3), 1)\n",
     "Ga0_3 = create_func(F,p_tot(1, 0.145), 1)\n",
     "\n",
     "c_tot_5 = create_func(F, p_tot(0.9, 0.18), 1)\n",
-    "c0_5 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.18), 1)\n",
+    "c0_5 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.18), 1)\n",
     "Ga0_5 = create_func(F,p_tot(1, 0.34), 1)\n",
     "\n",
     "c_tot_8 = create_func(F, p_tot(0.9, 0.1125), 1)\n",
-    "c0_8 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.1125), 1)\n",
+    "c0_8 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.1125), 1)\n",
     "Ga0_8 = create_func(F,p_tot(1, 0.8), 1)\n",
     "\n",
     "c_tot_9 = create_func(F, p_tot(0.9, 0.1), 1)\n",
-    "c0_9 = create_func(F, 'x[0]<0.1 ? 0 :'+p_tot(0.9, 0.1), 1)\n",
+    "c0_9 = create_func(F, 'x[0]<'+str(bp)+'? 0 :' +p_tot(0.9, 0.1), 1)\n",
     "Ga0_9 = create_func(F,p_tot(1, 1), 1)\n",
     "\n",
-    "sym = 0\n",
-    "c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, 0.01, 10, sym)\n",
-    "c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, 0.01, 10, sym)\n",
-    "c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, 0.01, 10, sym)\n",
-    "c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, 0.01, 10, sym)\n",
-    "c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, 0.01, 10, sym)\n",
-    "c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, 0.01, 10, sym)"
+    "c0_1 = calc_sim(c0_1, c_tot_1, Ga0_1, dt, 10, sym)\n",
+    "c0_2 = calc_sim(c0_2, c_tot_2, Ga0_2, dt, 10, sym)\n",
+    "c0_3 = calc_sim(c0_3, c_tot_3, Ga0_3, dt, 10, sym)\n",
+    "c0_5 = calc_sim(c0_5, c_tot_5, Ga0_5, dt, 10, sym)\n",
+    "c0_8 = calc_sim(c0_8, c_tot_8, Ga0_8, dt, 10, sym)\n",
+    "c0_9 = calc_sim(c0_9, c_tot_9, Ga0_9, dt, 10, sym)"
    ]
   },
   {
@@ -160,16 +163,17 @@
    "outputs": [],
    "source": [
     "# Plot outside, check invariance\n",
-    "x = np.linspace(0.1, 1, 1000)\n",
+    "x = np.linspace(bp, 1, 1000)\n",
     "y1 = eval_func(c0_1, x)\n",
     "y2 = eval_func(c0_2, x)\n",
     "y3 = eval_func(c0_3, x)\n",
     "y9 = eval_func(c0_9, x)\n",
-    "plt.plot((x-np.min(x)), y1)\n",
-    "plt.plot((x-np.min(x))/2, 2*y2)\n",
-    "plt.plot((x-np.min(x))/3, 3*y3)\n",
-    "plt.plot((x-np.min(x))/9, 9*y9)\n",
-    "# plt.xlim(0.0, 0.951252)\n",
+    "plt.plot((x-bp), y1)\n",
+    "plt.plot((x-bp)/2, 2*y2)\n",
+    "plt.plot((x-bp)/3, 3*y3)\n",
+    "plt.plot((x-bp)/9, 9*y9)\n",
+    "plt.xlim(0.0, 0.02)\n",
+    "plt.ylim(0.2, 1)\n",
     "plt.show()"
    ]
   },
@@ -185,11 +189,11 @@
     "y2 = eval_func(c0_2, x)-np.min(eval_func(c0_2, x))\n",
     "y3 = eval_func(c0_3, x)-np.min(eval_func(c0_3, x))\n",
     "y9 = eval_func(c0_9, x)-np.min(eval_func(c0_9, x))\n",
-    "plt.plot((x-np.min(x)), y1)\n",
-    "plt.plot((x-np.min(x)), y2)\n",
-    "plt.plot((x-np.min(x)), y3)\n",
-    "plt.plot((x-np.min(x)), y9)\n",
-    "plt.xlim(0.0, 0.1251252)\n",
+    "plt.plot((x-np.min(x)), y1/eval_func(c0_1, [0.08]))\n",
+    "plt.plot((x-np.min(x)), y2/eval_func(c0_2, [0.08]))\n",
+    "plt.plot((x-np.min(x)), y3/eval_func(c0_3, [0.08]))\n",
+    "plt.plot((x-np.min(x)), y9/eval_func(c0_9, [0.08]))\n",
+    "plt.xlim(0.08, bp)\n",
     "plt.show()"
    ]
   },