diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90
index d146046..2a4a331 100644
--- a/plasmaExpansion.f90
+++ b/plasmaExpansion.f90
@@ -307,6 +307,8 @@ program plasmaExpansion
 
   real(dp):: Temp_bc ! Temperature
   real(dp):: Temp0, TempF
+  real(dp):: Zave_bc ! Average charge state
+  real(dp):: Zave0, ZaveF
   real(dp):: n_ecr   ! Electron critical density for the laser
   real(dp):: c_s     ! Ion sound speed
   real(dp):: u_bc    ! Injection velocity
@@ -352,12 +354,14 @@ program plasmaExpansion
   ! Set input parameters (remember these have to be in non-dimensional units)
   Temp0    = 60.0_dp * eV_to_K / Temp_ref
   TempF    = 10.0_dp * eV_to_K / Temp_ref
+  Zave0    = 12.0_dp
+  ZaveF    =  3.0_dp
   n_ecr    = 1.0e19_dp * cm3_to_m3 / n_ref
-  c_s      = sqrt(T_to_Z(Temp0) * gam * Temp0)
+  c_s      = sqrt(Zave0 * gam * Temp0)
   u_bc0    = c_s!sqrt(Temp0)
   u_bcF    = sqrt(TempF)
-  n_bc0    = n_ecr / T_to_Z(Temp0)
-  n_bcF    = n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0)
+  n_bc0    = n_ecr / Zave0
+  n_bcF    = 1.0e-1*n_bc0!n_ecr*1.0e-1 / T_to_Z(Temp0)
 
   ! Set domain boundaries (non-dimensional units)
   r0 = 200.0e-6_dp / L_ref
@@ -487,19 +491,22 @@ program plasmaExpansion
 
   ! Main loop
   Temp_bc = Temp0
+  Zave_bc = Zave0
   u_bc    = u_bc0
   n_bc    = n_bc0
   do t = 1, nt
     if      (t >  t_bc0 .and. t <= t_bcF) then
       Temp_bc = (TempF - Temp0) / float(t_bcF - t_bc0)*(t - t_bc0) + Temp0
+      Zave_bc = (ZaveF - Zave0) / float(t_bcF - t_bc0)*(t - t_bc0) + Zave0
       u_bc    = (u_bcF - u_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + u_bc0
       n_bc    = (n_bcF - n_bc0) / float(t_bcF - t_bc0)*(t - t_bc0) + n_bc0
     else if (t >  t_bcF) then
       Temp_bc = TempF
+      Zave_bc = ZaveF
       u_bc    = u_bcF
       n_bc    = n_bcF
     end if
-    call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, T_to_Z(Temp_bc))
+    call writeOutputBoundary(t, dt, n_bc, u_bc, Temp_bc, Zave_bc)
     f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
     f0 = f0 * n_bc / (sum(f0)*dv)
 
@@ -508,7 +515,7 @@ program plasmaExpansion
     f_i_old(1,j0:nv) = f0
     f_i(1,j0:nv) = f_i_old(1,j0:nv)
     T_i(1) = Temp_bc
-    Zave(1) = T_to_Z(Temp_bc)
+    Zave(1) = Zave_bc
     ! r = rf, v<0
     f_i_old(nr,1:j0-1) = 0.0_dp
     f_i(nr,1:j0-1) = f_i_old(nr,1:j0-1)
@@ -534,7 +541,7 @@ program plasmaExpansion
         u_i(i)  = sum(v(:)   *f_i(i,:))*dv / n_i(i)
         E_i(i)  = sum(v(:)**2*f_i(i,:))*dv / n_i(i)
         T_i(i)  = 2.0_dp*E_i(i) - 2.0_dp*u_i(i)**2
-        Zave(i) = T_to_Z(T_i(i))
+        Zave(i) = Zave_bc
 
       else
         u_i(i)  = 0.0_dp
@@ -574,7 +581,7 @@ program plasmaExpansion
 
       ! Iterate system
       call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info)
-      phi = phi_old + Res
+      phi = phi_old + 1.0e-1_dp*Res
       ! phi0=phi(1)
 
       ! Calculate distribution of electrons