diff --git a/src/pyfmi/fmi.pyx b/src/pyfmi/fmi.pyx
index 5204cbd0..f4de45d8 100644
--- a/src/pyfmi/fmi.pyx
+++ b/src/pyfmi/fmi.pyx
@@ -8304,7 +8304,6 @@ 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)
@@ -8363,27 +8362,11 @@ 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]), 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]
+                eps_pt[i] = RUROUND*(max(abs(v_pt[i]), nominals_pt[i]))
         else:
             for i in range(len_v):
-                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]
-
+                tmp_nominal = self.get_variable_nominal(valueref = v_ref_pt[i])
+                eps_pt[i] = RUROUND*(max(abs(v_pt[i]), tmp_nominal))
 
         if group is not None:
             if output_matrix is not None: