diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex
index 75f0106..f9f9ae8 100644
--- a/doc/user-manual/fpakc_UserManual.tex
+++ b/doc/user-manual/fpakc_UserManual.tex
@@ -697,7 +697,8 @@ make
                                \begin{itemize}
                                \item \textbf{constant}: The flow does not change.
                                  This is the default value if none is provided.
-                               \item \textbf{quasiNeutral}: The weight of the injected species' particles is changed over time to keep a zero electric field at the boundary.
+                               \item \textbf{quasiNeutral}: The flow of the species is changed over time to keep the charge density close to zero at the injection surface.
+                                 The weight of the injected particles is changed.
                                \end{itemize}
           \item \textbf{physicalSurface}: Integer.
                                           Identification of the edge in the mesh file.
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 2b23ac6..942189c 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -53,7 +53,6 @@ MODULE moduleInject
   end type injectconstant
 
   type, extends(injectGeneric):: injectQuasiNeutral
-    real(8):: qSign_incident
 
   end type injectQuasiNeutral
 
@@ -229,14 +228,6 @@ MODULE moduleInject
       !Scale particles for different species steps
       IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
 
-      ! Specific parameters per type
-      select type(self)
-      type is(injectQuasiNeutral)
-        ! Get the sign of the charged species
-        self%qSign_incident = sign(1.d0, self%species%qm)
-
-      end select
-
     END SUBROUTINE initInject
 
     ! Update the value of injects
@@ -257,15 +248,17 @@ MODULE moduleInject
 
     ! Updates the flow in the injection to maintain quasineutrality
     subroutine updateQuasiNeutral(self)
-      use moduleMesh, only: meshEdge, meshCell
+      use moduleMesh, only: meshEdge, meshNode, mesh, qSpecies
       implicit none
 
       class(injectGeneric), intent(inout):: self
-      integer:: e
+      integer:: e, s, n
+      integer, allocatable:: nodes(:)
       class(meshEdge), pointer:: edge
-      class(meshCell), pointer:: cell
-      real(8):: Xi(1:3)
-      real(8):: EF_normal
+      real(8), allocatable:: density_nodes(:)
+      class(meshNode), pointer:: node
+      real(8):: den_center
+      real(8):: density_incident, density_rest
       real(8):: alpha
 
       select type(self)
@@ -273,21 +266,42 @@ MODULE moduleInject
         do e = 1, self%nEdges
           edge => self%edges(e)%obj
 
-          Xi = edge%centerXi()
+          density_incident = 0.d0
+          density_rest     = 0.d0
 
-          if (associated(edge%e1)) then
-            cell => edge%e1
+          nodes = edge%getNodes(edge%nNodes)
+          allocate(density_nodes(1:edge%nNodes))
+          do s = 1, nSpecies
+            do n = 1, edge%nNodes
+              node => mesh%nodes(nodes(n))%obj
 
-          else
-            cell => edge%e2
+              density_nodes(n) = node%output(s)%den
 
-          end if
+            end do
 
-          ! Projection of EF on the edge normal vector
-          EF_normal = dot_product(cell%gatherElectricField(Xi), edge%normal)
+            ! Gather the density at the edge center and multiply by the species charge
+            den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
+
+            if (s == self%species%n) then
+              density_incident = den_center
+
+            else
+              density_rest = density_rest + den_center
+
+            end if
+
+          end do
 
           ! Correction for this time step
-          alpha = self%qSign_incident * EF_normal
+          if (density_rest > 1.0d-10) then
+            ! If there is a rest population, correct
+            alpha = 1.d0 + density_incident/density_rest
+
+          else
+            ! If not, alpha is assumed unchaged
+            alpha = 0.d0
+
+          end if
 
           ! Adjust the weight of particles to match the new current
           self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-1 * alpha
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index 635c9eb..e0ca06c 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -466,7 +466,7 @@ MODULE moduleSolver
       use moduleInject
       IMPLICIT NONE
 
-      INTEGER:: i, sp
+      INTEGER:: i, n, sp
       
       DO i=1, nInject
         associate(inject => injects(i)%obj)