diff --git a/FloryHugg_DiffUnbleached.ipynb b/FloryHugg_DiffUnbleached.ipynb
index f9d2427fbc75ba1e75a53285192b6dca219d1883..d631829cc624d2d1bce084089a437a6e20aa0ed3 100644
--- a/FloryHugg_DiffUnbleached.ipynb
+++ b/FloryHugg_DiffUnbleached.ipynb
@@ -16,7 +16,7 @@
     "# domain = ms.Sphere(df.Point(0, 0, 0), 1.0)\n",
     "# mesh = ms.generate_mesh(domain, 50)\n",
     "mesh = df.UnitIntervalMesh(10000)\n",
-    "dt = 0.0001\n",
+    "dt = 0.001\n",
     "\n",
     "F = df.FunctionSpace(mesh, 'CG', 1)"
    ]
@@ -33,28 +33,28 @@
     "    X = df.SpatialCoordinate(mesh)\n",
     "#     X.interpolate(df.Expression('x[0]', degree=1))\n",
     "    # Weak form 1D:\n",
-    "#     form = ((df.inner((c-c0)/dt, tc) +\n",
-    "#              df.inner(df.grad(c), df.grad((1-c_tot)/Ga0*tc))) -\n",
-    "#              df.inner(df.grad(c_tot), df.grad((1-c_tot)/Ga0/c_tot*c*tc))-\n",
-    "#              tc*df.inner(df.grad(c), df.grad((1-c_tot)/Ga0))+\n",
-    "#              tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c/Ga0))) * df.dx\n",
+    "    form = ((df.inner((c-c0)/dt, tc) +\n",
+    "             df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc))) -\n",
+    "             df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc))-\n",
+    "             tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0))+\n",
+    "             tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0))) * df.dx\n",
     "    \n",
     "#     # Weak form radial symmetry:\n",
-    "    form = ((df.inner((c-c0)/dt, tc*X[0]*X[0]) +\n",
-    "             df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc*X[0]*X[0]))) -\n",
-    "             df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc*X[0]*X[0]))-\n",
-    "             tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*X[0]*X[0]))+\n",
-    "             tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0*X[0]*X[0]))-\n",
-    "             (1-c_tot)*Ga0*2*X[0]*c.dx(0)*tc+\n",
-    "             (1-c_tot)*Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx\n",
+    "#     form = ((df.inner((c-c0)/dt, tc*X[0]*X[0]) +\n",
+    "#              df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*tc*X[0]*X[0]))) -\n",
+    "#              df.inner(df.grad(c_tot), df.grad((1-c_tot)*Ga0/c_tot*c*tc*X[0]*X[0]))-\n",
+    "#              tc*df.inner(df.grad(c), df.grad((1-c_tot)*Ga0*X[0]*X[0]))+\n",
+    "#              tc*df.inner(df.grad(c_tot), df.grad((1-c_tot)/c_tot*c*Ga0*X[0]*X[0]))-\n",
+    "#              (1-c_tot)*Ga0*2*X[0]*c.dx(0)*tc+\n",
+    "#              (1-c_tot)*Ga0/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx\n",
     "  # Weak form radial symmetry:\n",
     "#     form = ((df.inner((c-c0)/dt, tc*X[0]*X[0]) +\n",
-    "#              df.inner(df.grad(c), df.grad((1-c_tot+a*c_tot*c_tot)*tc*X[0]*X[0]))) -\n",
-    "#              df.inner(df.grad(c_tot), df.grad((1-c_tot+a*c_tot*c_tot)/c_tot*c*tc*X[0]*X[0]))-\n",
-    "#              tc*df.inner(df.grad(c), df.grad((1-c_tot+a*c_tot*c_tot)*X[0]*X[0]))+\n",
-    "#              tc*df.inner(df.grad(c_tot), df.grad((1-c_tot+a*c_tot*c_tot)/c_tot*c*X[0]*X[0]))-\n",
-    "#              (1-c_tot+a*c_tot*c_tot)*2*X[0]*c.dx(0)*tc+\n",
-    "#              (1-c_tot+a*c_tot*c_tot)/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx\n",
+    "#              df.inner(df.grad(c), df.grad((1-c_tot+Ga0*c_tot*c_tot)*tc*X[0]*X[0]))) -\n",
+    "#              df.inner(df.grad(c_tot), df.grad((1-c_tot+Ga0*c_tot*c_tot)/c_tot*c*tc*X[0]*X[0]))-\n",
+    "#              tc*df.inner(df.grad(c), df.grad((1-c_tot+Ga0*c_tot*c_tot)*X[0]*X[0]))+\n",
+    "#              tc*df.inner(df.grad(c_tot), df.grad((1-c_tot+Ga0*c_tot*c_tot)/c_tot*c*X[0]*X[0]))-\n",
+    "#              (1-c_tot+Ga0*c_tot*c_tot)*2*X[0]*c.dx(0)*tc+\n",
+    "#              (1-c_tot+Ga0*c_tot*c_tot)/c_tot*c*2*X[0]*c_tot.dx(0)*tc) * df.dx\n",
     "    \n",
     "    t = 0\n",
     "    # Solve in time\n",
@@ -71,7 +71,11 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    }
+   },
    "outputs": [],
    "source": [
     "# Interpolate c_tot and initial conditions\n",
@@ -110,6 +114,42 @@
     "Ga0_9.interpolate(df.Expression('1', degree=1))"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p1_i = 0.9\n",
+    "p1_o = 0.1\n",
+    "p2_i = 0.8\n",
+    "p2_o = 0.2\n",
+    "a1 = 0\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",
+    "\n",
+    "ct_1 = df.Function(F)\n",
+    "ct_2 = df.Function(F)\n",
+    "c0_1 = df.Function(F)\n",
+    "c0_2 = df.Function(F)\n",
+    "g_1 = df.Function(F)\n",
+    "g_2 = df.Function(F)\n",
+    "\n",
+    "ct_1.interpolate(df.Expression(p_tot(p1_i, p1_o), degree=1))\n",
+    "ct_2.interpolate(df.Expression(p_tot(p2_i, p2_o), degree=1))\n",
+    "P1 = c_tot1(0)/c_tot1(1)\n",
+    "P2 = c_tot2(0)/c_tot2(1)\n",
+    "\n",
+    "g_1.interpolate(df.Expression('1', degree=1))\n",
+    "g_2.interpolate(df.Expression(p_tot(p1_o/p2_o, 1), degree=1))\n",
+    "D_out1 = 1-c_tot1(1)\n",
+    "a = (P2*D_out1/P1-1+p2_o)/p2_o**2\n",
+    "\n",
+    "c0_1.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p1_i, p1_o), degree=1))\n",
+    "c0_2.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p2_i, p2_o), degree=1))"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -133,9 +173,9 @@
    "outputs": [],
    "source": [
     "# 1D:\n",
-    "plt.plot(np.linspace(0, 1, 10000), [1.31*c0_1([x]) for x in np.linspace(0, 1, 10000)])\n",
+    "plt.plot(np.linspace(0, 1, 10000), [1.465*c0_1([x]) for x in np.linspace(0, 1, 10000)])\n",
     "plt.plot(np.linspace(0, 1, 10000), [c0_2([x]) for x in np.linspace(0, 1, 10000)])\n",
-    "plt.xlim(0.06, 0.125)\n",
+    "plt.xlim(0.0, 0.125)\n",
     "# plt.ylim(0, 0.3)\n",
     "# 3D:\n",
     "# plt.plot(np.linspace(0, 0.5, 1000), [c0([x, 0, 0]) for x in np.linspace(0, 0.5, 1000)])"
@@ -166,44 +206,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "p1_i = 0.9\n",
-    "p1_o = 0.1\n",
-    "p2_i = 0.8\n",
-    "p2_o = 0.2\n",
-    "a1 = 0\n",
-    "\n",
-    "def p_tot(p_i, p_o):\n",
-    "    return str(p_i-p_o)+'*(-0.5*tanh(350000*(x[0]-0.1))+0.5)+'+str(p_o)\n",
-    "\n",
-    "ct_1 = df.Function(F)\n",
-    "ct_2 = df.Function(F)\n",
-    "c0_1 = df.Function(F)\n",
-    "c0_2 = df.Function(F)\n",
-    "g_1 = df.Function(F)\n",
-    "g_2 = df.Function(F)\n",
-    "\n",
-    "ct_1.interpolate(df.Expression(p_tot(p1_i, p1_o), degree=1))\n",
-    "ct_2.interpolate(df.Expression(p_tot(p2_i, p2_o), degree=1))\n",
-    "P1 = c_tot1(0)/c_tot1(1)\n",
-    "P2 = c_tot2(0)/c_tot2(1)\n",
-    "\n",
-    "g_1.interpolate(df.Expression('1', degree=1))\n",
-    "g_2.interpolate(df.Expression(p_tot(p1_o/p2_o, p1_i/p2_i), degree=1))\n",
-    "D_out1 = 1-c_tot1(1)\n",
-    "a = (P2*D_out1/P1-1+p2_o)/p2_o**2\n",
-    "\n",
-    "c0_1.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p1_i, p1_o), degree=1))\n",
-    "c0_2.interpolate(df.Expression('x[0]<0.1 ? 0 :'+p_tot(p2_i, p2_o), degree=1))"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.plot(np.linspace(0, 1, 1000), [ct_1([x]) for x in np.linspace(0, 1, 1000)])\n",
-    "plt.plot(np.linspace(0, 1, 1000), [ct_2([x]) for x in np.linspace(0, 1, 1000)])\n",
+    "plt.plot(np.linspace(0, 1, 1000), [(1-ct_1([x]))*ct_1([x]) for x in np.linspace(0, 1, 1000)])\n",
+    "plt.plot(np.linspace(0, 1, 1000), [(1-ct_2([x]))*g_2([x]) for x in np.linspace(0, 1, 1000)])\n",
     "plt.show()\n",
     "\n",
     "plt.plot(np.linspace(0, 1, 1000), [g_2([x]) for x in np.linspace(0, 1, 1000)])\n",