diff --git a/plasmaExpansion.f90 b/plasmaExpansion.f90
index f715510..896629d 100644
--- a/plasmaExpansion.f90
+++ b/plasmaExpansion.f90
@@ -35,7 +35,7 @@ program plasmaExpansion
   real(dp), parameter:: m_i     = 1.9712258e-25_dp ! Tin atom mass in kg
   real(dp), parameter:: gamma_i = 1.0_dp           ! Adiabatic coefficient for ions
   real(dp), parameter:: m_e     = 9.1093837e-31_dp ! Electron mass in kg
-  real(dp), parameter:: gamma_e = 2.0_dp           ! Adiabatic coefficient for electrons
+  real(dp), parameter:: gamma_e = 4.0_dp / 3.0_dp  ! Adiabatic coefficient for electrons
 
   real(dp):: r0, rf
   real(dp), allocatable, dimension(:):: r
@@ -285,9 +285,9 @@ program plasmaExpansion
       phi_old = phi
 
       ! Diagonal matrix for Newton integration scheme
-      db_dphi        = n_e / T_e ! Isotropic
-      ! db_dphi        =  Zave(1) * n_i(1) / (gamma_e * T_e) * &
-                      ! (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi_old-phi0)/T_e)**((2.0_dp - gamma_e)/(gamma_e - 1.0_dp))) ! Polytropic
+      ! db_dphi        = n_e / T_e ! Isotropic
+      db_dphi        =  Zave(1) * n_i(1) / (gamma_e * T_e) * &
+                      (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi_old-phi0)/T_e)**((2.0_dp - gamma_e)/(gamma_e - 1.0_dp)) ! Polytropic
       diag           = -2.0_dp / dr**2 - db_dphi
       diag_low       =  1.0_dp / dr**2 - 1.0_dp / (r(2:nr)   * dr)
       diag_high      =  1.0_dp / dr**2 + 1.0_dp / (r(1:nr-1) * dr)
@@ -314,8 +314,8 @@ program plasmaExpansion
       ! phi0=phi(1) ! Neumann
 
       ! Calculate distribution of electrons
-      n_e = Zave(1) * n_i(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
-      ! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Polytropic
+      ! n_e = Zave(1) * n_i(1) * exp((phi- phi0) / T_e) ! Isothermal (Boltzmann)
+      n_e = Zave(1) * n_i(1) * (1.0_dp + (gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp)) ! Polytropic
 
       ! Check if the solution has converged
       phiConv = maxval(abs(Res),1)