diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index 492eaf5..e9742c9 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -1306,6 +1306,9 @@ MODULE moduleInput
 
           end do
 
+          allocate(bound%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
+          bound%alpha = 0.d0
+
         end select
 
       end do
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 25d7175..1cb4dc7 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -651,6 +651,7 @@ MODULE moduleMesh
     integer:: n
     character(:), allocatable:: name 
     procedure(updateParticle_interface), pointer, pass:: update => null()
+    procedure(printParticle_interface),  pointer, pass:: print  => null()
     contains
       procedure(applyParticle_interface), deferred, pass:: apply
 
@@ -708,6 +709,14 @@ MODULE moduleMesh
 
     end subroutine updateParticle_interface
 
+    ! Update the values of the particle boundary model
+    subroutine printParticle_interface(self)
+      import boundaryParticleGeneric
+
+      class(boundaryParticleGeneric), intent(inout):: self
+
+    end subroutine printParticle_interface
+
   end interface
 
   !Reflecting boundary
diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90
index c51c26d..c30ca1d 100644
--- a/src/modules/mesh/moduleMesh@boundaryParticle.f90
+++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90
@@ -109,10 +109,9 @@ submodule(moduleMesh) boundaryParticle
 
         boundary%s_incident = s_incident
 
-        allocate(boundary%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
-        boundary%alpha = 0.d0
 
         boundary%update => quasiNeutrality_update
+        boundary%print  => quasiNeutrality_print
 
       end select
 
@@ -391,7 +390,7 @@ submodule(moduleMesh) boundaryParticle
           alpha = 1.d0 - density_incident/density_rest
 
           ! Apply correction with a factor of 0.1 to avoid fast changes
-          self%alpha(edge%n) = self%alpha(edge%n) + 0.1d0 * alpha
+          self%alpha(edge%n) = self%alpha(edge%n) + 1.0d-2 * alpha
 
           ! Limit alpha between 0 and 1
           self%alpha(edge%n) = min(self%alpha(edge%n), 1.d0)
@@ -405,6 +404,21 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine quasiNeutrality_update
 
+    subroutine quasiNeutrality_print(self)
+      implicit none
+
+      class(boundaryParticleGeneric), intent(inout):: self
+
+      print *, 'test'
+
+      select type(self)
+      type is(boundaryQuasiNeutrality)
+        print*, self%alpha
+
+      end select
+
+    end subroutine quasiNeutrality_print
+
     ! Generic boundary conditions for internal use
     module subroutine genericReflection(edge, part)
       use moduleCaseParam
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index cfe473a..4ea051c 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -523,6 +523,8 @@ MODULE moduleSolver
       USE moduleCaseParam, ONLY: timeStep
       IMPLICIT NONE
 
+      integer:: b ! TEMPORARY
+
       CALL outputProbes()
 
       counterOutput = counterOutput + 1
@@ -551,6 +553,17 @@ MODULE moduleSolver
         END IF
         WRITE(*,*)
 
+        ! TEMPORARY
+        ! Output of boundaries. TODO: Move this to mesh module
+        do b = 1, nBoundariesParticle
+          if (associated(boundariesParticle(b)%obj%print)) then
+            print *, b
+            call boundariesParticle(b)%obj%print()
+
+          end if
+
+        end do
+
       END IF
 
       counterCPUTime = counterCPUTime + 1