diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index da43e98..e68a732 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -314,6 +314,7 @@ MODULE moduleInput
       USE moduleOutput
       USE moduleRefParam
       USE moduleRandom
+      USE moduleProbe
       USE json_module
       IMPLICIT NONE
 
@@ -420,6 +421,9 @@ MODULE moduleInput
 
               END IF
 
+              !Assign particles to probes
+              CALL doProbes(partNew)
+
             END DO
 
             DEALLOCATE(sourceScalar, sourceArray)
@@ -1089,6 +1093,8 @@ MODULE moduleInput
 
       END DO
 
+      CALL resetProbes(tInitial)
+
     END SUBROUTINE readProbes
 
     SUBROUTINE readEMBoundary(config)
diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
index d50263b..4341cb0 100644
--- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90
+++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
@@ -508,7 +508,7 @@ MODULE moduleMesh2DCart
       conv = 1.D0
       XiO  = 0.D0
 
-      DO WHILE(conv > 1.D-2)
+      DO WHILE(conv > 1.D-4)
         dPsi = self%dPsi(XiO, 4)
         pDer = self%partialDer(4, dPsi)
         detJ = self%detJac(pDer)
@@ -516,8 +516,7 @@ MODULE moduleMesh2DCart
         fPsi = self%fPsi(XiO, 4)
         f = (/ DOT_PRODUCT(fPsi,self%x), &
                DOT_PRODUCT(fPsi,self%y), &
-               0.D0 /)
-        f = f - r
+               0.D0 /) - r
         Xi = XiO - MATMUL(invJ, f)/detJ
         conv = MAXVAL(DABS(Xi-XiO),1)
         XiO = Xi
diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
index f8e41b3..a2ffb7a 100644
--- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
@@ -532,8 +532,7 @@ MODULE moduleMesh2DCyl
         fPsi = self%fPsi(XiO, 4)
         f = (/ DOT_PRODUCT(fPsi,self%z), &
                DOT_PRODUCT(fPsi,self%r), &
-               0.D0 /)
-        f = f - r
+               0.D0 /) - r
         Xi = XiO - MATMUL(invJ, f)/detJ
         conv = MAXVAL(DABS(Xi-XiO),1)
         XiO = Xi
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index cb252e8..9a3c4ef 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -872,7 +872,7 @@ MODULE moduleMesh
 
                   !If particles do not belong to the species, skip collision
                   !This can happen, for example, if particle has been previously ionized or removed
-                  !TODO: Try to find a way to no lose these collisions. Maybe check new 'k' and use that for the collision, maybe?
+                  !TODO: Try to find a way to not lose these collisions. Maybe check new 'k' and use that for the collision?
                   IF (part_i%species%n /= i .OR. &
                       part_j%species%n /= j) THEN
                       CYCLE
diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90
index 2d2295a..c7d3cf5 100644
--- a/src/modules/moduleProbe.f90
+++ b/src/modules/moduleProbe.f90
@@ -91,11 +91,17 @@ MODULE moduleProbe
                       1:self%nv(3)))
 
       !Number of iterations between output
-      self%every = NINT(timeStep/ tauMin / ti_ref)
+      IF (timeStep == 0.D0) THEN
+        self%every = 1
+
+      ELSE
+        self%every = NINT(timeStep/ tauMin / ti_ref)
+
+      END IF
       
       !Maximum radius
       !TODO: Make this an input parameter
-      self%maxR = 1.D-1/L_ref
+      self%maxR = 1.D0
 
       !Init the probe lock
       CALL OMP_INIT_LOCK(self%lock)
@@ -110,6 +116,8 @@ MODULE moduleProbe
       INTEGER, INTENT(out):: i, j, k
       LOGICAL, INTENT(out):: inside
 
+      inside = .TRUE.
+
       i = FLOOR((vp(1) - self%vi(1))/self%vrange(1)*(REAL(self%nv(1) - 1)) + 1.D0)
       IF (i >= self%nv(1) .OR. i <  1) inside = .FALSE.
       j = FLOOR((vp(2) - self%vj(1))/self%vrange(2)*(REAL(self%nv(2) - 1)) + 1.D0)
@@ -198,7 +206,7 @@ MODULE moduleProbe
 
       WRITE(tstring, iterationFormat) t
       WRITE(pstring, "(I3.3)") self%id
-      fileName='OUTPUT_' // tstring// '_f_' // pstring // '.dat'
+      fileName='Probe_' // tstring// '_f_' // pstring // '.dat'
       WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
       OPEN (10, file = path // folder // '/' //  fileName)
       WRITE(10, "(A1, 1X, A)")               "# ", self%species%name
@@ -268,7 +276,7 @@ MODULE moduleProbe
 
       DO i = 1, nProbes
         probe(i)%f = 0.D0
-        probe(i)%update = MOD(t, probe(i)%every) == 0 .OR. t == tFinal
+        probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0
 
       END DO
 
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index da2c9c4..69d4055 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -523,6 +523,8 @@ MODULE moduleSolver
 
       END IF
 
+      CALL outputProbes(t)
+
       counterOutput = counterOutput + 1
       IF (counterOutput >= triggerOutput .OR. &
           t == tFinal .OR. t == tInitial) THEN
@@ -530,7 +532,6 @@ MODULE moduleSolver
         !Resets output counter
         counterOutput=0
 
-        CALL outputProbes(t)
         CALL mesh%printOutput(t)
         IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
         CALL mesh%printEM(t)