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 83a567a..c7d3cf5 100644
--- a/src/modules/moduleProbe.f90
+++ b/src/modules/moduleProbe.f90
@@ -91,7 +91,13 @@ 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
@@ -270,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 .OR. t == tInitial
+        probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0
 
       END DO