diff --git a/bc_5ns.csv b/bc_5ns.csv
index f364754..306b9b5 100644
--- a/bc_5ns.csv
+++ b/bc_5ns.csv
@@ -1,6 +1,5 @@
 t (s),n (m^-3),u (m s^-1),T (eV),Z
 0.000000E-000,1.111111E+025,0.000000E+000,3.000000E+001,9.000000E+000
-1.000000E-008,1.111111E+025,0.000000E+000,3.000000E+001,9.000000E+000
-1.500000E-008,1.000000E+022,0.000000E+000,5.000000E+000,9.000000E+000
-2.000000E-008,1.000000E+016,0.000000E+000,5.000000E+000,9.000000E+000
-2.000000E-006,1.000000E+010,0.000000E+000,5.000000E+000,9.000000E+000
+5.000000E-008,1.111111E+025,0.000000E+000,3.000000E+001,9.000000E+000
+5.500000E-008,1.000000E+020,0.000000E+000,5.000000E+000,9.000000E+000
+2.000000E-006,1.000000E+020,0.000000E+000,5.000000E+000,9.000000E+000
diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90
index 20b456c..fd24795 100644
--- a/plasmaExpansion.f90
+++ b/plasmaExpansion.f90
@@ -100,7 +100,7 @@ program plasmaExpansion
   ! Set domain boundaries (non-dimensional units)
   r0 =  10.0e-6_dp / L_ref
   rf =   2.0e-3_dp / L_ref
-  dr = 5.0e2_dp
+  dr =   1.0e-6_dp / L_ref
   nr = nint((rf - r0) / dr) + 1
   dr =      (rf - r0) / float(nr-1)
   allocate(r(1:nr))
@@ -280,7 +280,7 @@ program plasmaExpansion
     n_e = Zave * n_i
 
     ! Solve Poission (maximum number of iterations, break if convergence is reached before)
-    do k = 1, 1000
+    do k = 1, 2000
       ! Store previous value
       phi_old = phi
 
@@ -319,7 +319,7 @@ program plasmaExpansion
 
       ! Check if the solution has converged
       phiConv = maxval(abs(Res),1)
-      if (phiConv < 1.0e-4_dp) then
+      if (phiConv < 1.0e-6_dp) then
         exit
 
       end if
diff --git a/scripts_python/plotCumF.py b/scripts_python/plotCumF.py
index 189ef2a..e0896b5 100644
--- a/scripts_python/plotCumF.py
+++ b/scripts_python/plotCumF.py
@@ -11,27 +11,28 @@ m_i = 1.9712e-25
 #  paths = ['../quasiNeutral_partialAblation/','../Poisson_partialAblation/','../PoissonTi_partialAblation/','../2024-11-05_13.56.23/']
 #  paths = ['../2024-11-11_13.58.54/']
 #  paths = ['../quasiNeutral_fullAblation/','../Poisson_fullAblation/', '../quasiNeutral_partialAblation/', '../Poisson_partialAblation/']
-paths = ['../2024-11-18_11.54.30/']
+paths = ['../2024-11-28_10.01.56/']
+#  paths = ['../Poisson_partialAblation/','../Poisson_partialAblation_lowerT/','../Poisson_partialAblation_lowT/','../Poisson_partialAblation_highT/']
 labels = [path[3:-1] for path in paths]
 
 for path, label in zip(paths, labels):
     filesCum_i = sorted(glob.glob(path+'time_*_fCum_i.csv'))
     start = 0
     end   = len(filesCum_i)
-    every = 50
+    every = 20
     for fileCum_i in filesCum_i[start:end+1:every]:
         time, x, v, f_i = readF.read(fileCum_i)
         plt.plot(v**2*m_i*0.5/e, f_i[0]*e/m_i/v, label='t = {:.1f} ns'.format(time*1e9))
 
-    #  time, x, v, f_i = readF.read(filesCum_i[-1])
-    #  plt.plot(v**2*m_i*0.5/e, f_i[0]*e/m_i/v, label=label)
+    time, x, v, f_i = readF.read(filesCum_i[-1])
+    plt.plot(v**2*m_i*0.5/e, f_i[0]*e/m_i/v, label=label)
 
-    plt.yscale('log')
-    plt.ylim([1e16,5e27])
-    plt.ylabel('Sum f(e) / sqrt(e) (m^-3 eV^-1)')
-    plt.xscale('log')
-    plt.xlim([1e0,1e4])
-    plt.xlabel('e (eV)')
+plt.yscale('log')
+plt.ylim([1e16,5e27])
+plt.ylabel('Sum f(e) / sqrt(e) (m^-3 eV^-1)')
+plt.xscale('log')
+plt.xlim([1e0,1e4])
+plt.xlabel('e (eV)')
 
 plt.legend()
 
diff --git a/scripts_python/plotMom.py b/scripts_python/plotMom.py
index 57e7dce..95ea48a 100644
--- a/scripts_python/plotMom.py
+++ b/scripts_python/plotMom.py
@@ -11,25 +11,25 @@ from scipy.constants import e, k
 #  paths = ['../quasiNeutral_partialAblation/','../Poisson_partialAblation/']
 #  paths = ['../2024-10-02_14.30.44/']
 #  paths = ['../quasiNeutral_fullAblation/','../Poisson_fullAblation/']
-paths = ['../2024-11-18_08.57.04/']
+paths = ['../2024-11-28_10.01.56/']
 labels = [path[3:-1] for path in paths]
 
 for path, label in zip(paths, labels):
     filesPhi   = sorted(glob.glob(path+'time_*_phi.csv'))
     filesMom_i = sorted(glob.glob(path+'time_*_mom_i.csv'))
-    start = 500
-    end   = len(filesMom_i)
-    every = 200
+    start = 0
+    end   = 50#len(filesMom_i)
+    every = 10
     fig, ax = plt.subplots(4, sharex='all')
     for fileMom_i, filePhi in zip(filesMom_i[start:end+1:every], filesPhi[start:end+1:every]):
         time, r, phi, E, n_e         = readPhi.read(filePhi)
         time, r, n_i, u_i, T_i, Zave = readMom.read(fileMom_i)
 
         #  ax[0].plot(r, phi, label='t = {:.1f} ns'.format(time*1e9))
-        ax[0].plot(r, E, label='t = {:.1f} ns'.format(time*1e9))
-        #  ax[0].plot(r, (Zave*n_i - n_e)/n_e, label='t = {:.1f} ns'.format(time*1e9))
+        #  ax[0].plot(r, E, label='t = {:.1f} ns'.format(time*1e9))
+        ax[0].plot(r, (Zave*n_i - n_e), label='t = {:.1f} ns'.format(time*1e9))
         ax[1].set_yscale('log')
-        ax[1].set_ylim([1e20,2e25])
+        #  ax[1].set_ylim(bottom=1e16)
         ax[1].plot(r, Zave*n_i)
         ax[1].plot(r,      n_e, color='k', linestyle='dashed')
         ax[2].plot(r,      u_i)