diff --git a/src/pyfmi/fmi.pyx b/src/pyfmi/fmi.pyx
index f4de45d8..5204cbd0 100644
--- a/src/pyfmi/fmi.pyx
+++ b/src/pyfmi/fmi.pyx
@@ -8304,6 +8304,7 @@ cdef class FMUModelME2(FMUModelBase2):
         cdef double tmp_nominal, fac, tmp
         cdef int method = FORWARD_DIFFERENCE if self.force_finite_differences is True or self.force_finite_differences == 0 else CENTRAL_DIFFERENCE
         cdef double RUROUND = FORWARD_DIFFERENCE_EPS if method == FORWARD_DIFFERENCE else CENTRAL_DIFFERENCE_EPS
+        cdef double RRUROUND = (RUROUND)**(0.5)
         cdef N.ndarray[FMIL.fmi2_real_t, ndim=1, mode='c'] dfpert, df, eps, nominals
         cdef N.ndarray[FMIL.fmi2_value_reference_t, ndim=1, mode='c'] v_ref = N.array(var_ref, copy=False, dtype = N.uint32)
         cdef N.ndarray[FMIL.fmi2_value_reference_t, ndim=1, mode='c'] z_ref = N.array(func_ref, copy=False, dtype = N.uint32)
@@ -8362,11 +8363,27 @@ cdef class FMUModelME2(FMUModelBase2):
                         nominals_pt[i] = self.get_variable_nominal(valueref = v_ref_pt[i])
 
             for i in range(len_v):
-                eps_pt[i] = RUROUND*(max(abs(v_pt[i]), nominals_pt[i]))
+                eps_pt[i] = RUROUND*(max(abs(v_pt[i]), RRUROUND))
+                if N.sign(v_pt[i]):
+                    eps_pt[i] *= N.sign(v_pt[i])                     
+                    if abs(v_pt[i])/2 <= abs(v_pt[i] + eps_pt[i]) <= 2*abs(v_pt[i]):
+                        temp = v_pt[i] + eps_pt[i]      
+                        eps_pt[i] = temp - v_pt[i]
+                    elif abs(v_pt[i])/2 <= abs(v_pt[i] - eps_pt[i]) <= 2*abs(v_pt[i]):
+                        temp = v_pt[i] - eps_pt[i]
+                        eps_pt[i] = temp - v_pt[i]
         else:
             for i in range(len_v):
-                tmp_nominal = self.get_variable_nominal(valueref = v_ref_pt[i])
-                eps_pt[i] = RUROUND*(max(abs(v_pt[i]), tmp_nominal))
+                eps_pt[i] = RUROUND*(max(abs(v_pt[i]), RRUROUND))
+                if N.sign(v_pt[i]):
+                    eps_pt[i] *= N.sign(v_pt[i])                     
+                    if abs(v_pt[i])/2 <= abs(v_pt[i] + eps_pt[i]) <= 2*abs(v_pt[i]):
+                        temp = v_pt[i] + eps_pt[i]      
+                        eps_pt[i] = temp - v_pt[i]
+                    elif abs(v_pt[i])/2 <= abs(v_pt[i] - eps_pt[i]) <= 2*abs(v_pt[i]):
+                        temp = v_pt[i] - eps_pt[i]
+                        eps_pt[i] = temp - v_pt[i]
+
 
         if group is not None:
             if output_matrix is not None: