Change in parameters, improve stability

Finally I've fixed all issues with the electric potential. I think these
conditions are good to show the Diko's peak.

plasmaExpansion.f90

@@ -99,8 +99,8 @@ program plasmaExpansion
 
   ! Set domain boundaries (non-dimensional units)
   r0 =  10.0e-6_dp / L_ref
-  rf = 100.0e-6_dp / L_ref
+  rf =   1.0e-3_dp / L_ref
-  dr = 2.0e1_dp
+  dr = 5.0e2_dp
   nr = nint((rf - r0) / dr) + 1
   dr =      (rf - r0) / float(nr-1)
   allocate(r(1:nr))
@@ -109,7 +109,7 @@ program plasmaExpansion
   end do
 
   ! Set position to calculate cumulative sum of f (non-dimensional units)
-  rCum = 80.0e-6 / L_ref
+  rCum =  5.0e-4 / L_ref
 
   ! Index for cumulative sum
   rCum_index = minloc(abs(r - rCum), 1)
@@ -131,7 +131,7 @@ program plasmaExpansion
   end if
 
   t0 = 0.0_dp
-  tf = 1.0e-7_dp / t_ref
+  tf = 5.0e-7_dp / t_ref
   ! tf = 1.0e1_dp * (rf - r0) / c_s
   dt = 1.0e-2_dp*dr/c_s
   nt = nint((tf - t0) / dt) 
@@ -308,7 +308,7 @@ program plasmaExpansion
 
       ! Iterate system
       call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info)
-      phi = phi_old + 1.0e-1_dp*Res
+      phi = phi_old + Res
       ! phi0=phi(1) ! Neumann
 
       ! Calculate distribution of electrons
@@ -317,7 +317,7 @@ program plasmaExpansion
 
       ! Check if the solution has converged
       phiConv = maxval(abs(Res),1)
-      if (phiConv < 1.0e-3_dp) then
+      if (phiConv < 1.0e-4_dp) then
         exit
 
       end if