diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90
index 091e52e..e81e796 100644
--- a/src/modules/mesh/moduleMeshBoundary.f90
+++ b/src/modules/mesh/moduleMeshBoundary.f90
@@ -77,10 +77,23 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE transparent
 
+    !Symmetry axis. Reflects particles.
+    !Although this function should never be called, it is set as a reflective boundary
+    !to properly deal with possible particles reaching a corner and selecting this boundary.
+    SUBROUTINE symmetryAxis(edge, part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+
+      CALL reflection(edge, part)
+
+    END SUBROUTINE symmetryAxis
+
     !Wall with temperature
     SUBROUTINE wallTemperature(edge, part)
       USE moduleSpecies
-      USE moduleBoundary
       USE moduleRandom
       IMPLICIT NONE
 
@@ -204,19 +217,28 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE ionization
 
-    !Symmetry axis. Reflects particles.
-    !Although this function should never be called, it is set as a reflective boundary
-    !to properly deal with possible particles reaching a corner and selecting this boundary.
-    SUBROUTINE symmetryAxis(edge, part)
-      USE moduleSpecies
-      IMPLICIT NONE
+    subroutine quasiNeutrality(edge, part)
+      use moduleRandom
+      implicit none
 
-      CLASS(meshEdge), INTENT(inout):: edge
-      CLASS(particle), INTENT(inout):: part
+      class(meshEdge), intent(inout):: edge
+      class(particle), intent(inout):: part
 
-      CALL reflection(edge, part)
+      select type(bound => edge%boundary%bTypes(part%species%n)%obj)
+      type is(boundaryQuasiNeutrality)
+        bound%alpha = 1.0d-1
 
-    END SUBROUTINE symmetryAxis
+        if (random() < bound%alpha) then
+          call reflection(edge, part)
+
+        else
+          call transparent(edge, part)
+
+        end if
+
+      end select
+
+    end subroutine quasiNeutrality
 
     !Points the boundary function to specific type
     SUBROUTINE pointBoundaryFunction(edge, s)
@@ -236,17 +258,20 @@ MODULE moduleMeshBoundary
       TYPE IS(boundaryTransparent)
         edge%fBoundary(s)%apply => transparent
 
+      TYPE IS(boundaryAxis)
+        edge%fBoundary(s)%apply => symmetryAxis
+
       TYPE IS(boundaryWallTemperature)
         edge%fBoundary(s)%apply => wallTemperature
 
       TYPE IS(boundaryIonization)
         edge%fBoundary(s)%apply => ionization
 
-      TYPE IS(boundaryAxis)
-        edge%fBoundary(s)%apply => symmetryAxis
+      type is(boundaryQuasiNeutrality)
+        edge%fBoundary(s)%apply => quasiNeutrality
 
       CLASS DEFAULT
-        CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')
+        CALL criticalError("Boundary type not defined", 'pointBoundaryFunction')
 
       END SELECT
 
diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90
index 0b76105..2b81491 100644
--- a/src/modules/moduleBoundary.f90
+++ b/src/modules/moduleBoundary.f90
@@ -26,6 +26,12 @@ MODULE moduleBoundary
 
   END TYPE boundaryTransparent
 
+  !Symmetry axis
+  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
+    CONTAINS
+
+  END TYPE boundaryAxis
+
   !Wall Temperature boundary
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature
     !Thermal velocity of the wall: square root(Wall temperature X specific heat)
@@ -47,12 +53,11 @@ MODULE moduleBoundary
 
   END TYPE boundaryIonization
 
-  !Symmetry axis
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
-    CONTAINS
-
-  END TYPE boundaryAxis
+  ! Ensures quasi-neutrality by changing the reflection coefficient
+  type, public, extends(boundaryGeneric):: boundaryQuasiNeutrality
+    real(8):: alpha ! Reflection parameter
 
+  end type boundaryQuasiNeutrality
   !Wrapper for boundary types (one per species)
   TYPE:: bTypesCont
     CLASS(boundaryGeneric), ALLOCATABLE:: obj
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 528cb91..3962986 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -257,7 +257,7 @@ MODULE moduleInject
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
       REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistMaxwellian(vTh = DSQRT(temperature/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
 
     END SUBROUTINE initVelDistMaxwellian
 
@@ -267,7 +267,7 @@ MODULE moduleInject
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
       REAL(8), INTENT(in):: temperature, m
 
-      velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m))
+      velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
 
     END SUBROUTINE initVelDistHalfMaxwellian
 
@@ -289,7 +289,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      v = sqrt(2.0)*self%vTh*randomMaxwellian()
+      v = self%vTh*randomMaxwellian()/sqrt(2.d0)
 
     END FUNCTION randomVelMaxwellian
 
@@ -302,7 +302,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      v = sqrt(2.0)*self%vTh*randomHalfMaxwellian()
+      v = self%vTh*randomHalfMaxwellian()
 
     END FUNCTION randomVelHalfMaxwellian
 
@@ -324,6 +324,7 @@ MODULE moduleInject
       USE moduleMesh
       USE moduleRandom
       USE moduleErrors
+      use moduleMeshBoundary, only: reflection
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
@@ -387,16 +388,20 @@ MODULE moduleInject
 
           partInj(n)%v = 0.D0
 
-          partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
-                                                   self%v(2)%obj%randomVel(), &
-                                                   self%v(3)%obj%randomVel() /)
+          do while(dot_product(partInj(n)%v, direction) <= 0.d0)
+            partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
+                                                    self%v(2)%obj%randomVel(), &
+                                                    self%v(3)%obj%randomVel() /)
 
-          !If injecting a no-drift distribution and velocity is negative, reflect
-          if ((self%vMod == 0.D0) .and. &
-              (dot_product(direction, partInj(n)%v) < 0.D0)) then
-            partInj(n)%v = - partInj(n)%v
+            !If while injecting a no-drift distribution the velocity is negative, reflect
+            if ((self%vMod == 0.D0) .and. &
+                (dot_product(direction, partInj(n)%v) < 0.D0)) then
+              call reflection(randomEdge, partInj(n))
+
+            end if
+
+          end do
 
-          end if
 
           !Obtain natural coordinates of particle in cell
           partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)