Use the default Method->EulerSum when expr is not analytic in the neighborhood of x₀:

Check:

An incorrect answer is obtained with Method->NIntegrate:

Here is a derivative where the default method works poorly:

The correct answer is:

In this case the expression is analytic, so Method->NIntegrate will work well:

Use Scale->s to capture the region of variation:

The scale of variation is around .01:

A value of Scale->s that is too large can be compensated by increasing the number of terms:

Use Scale to specify directional derivatives. The left and right derivatives of the nonanalytic function :

Check:

Complex directions may also be specified:

Check:

Use the option Scale to avoid regions of non-analyticity when the method used is NIntegrate:

Shrinking the radius avoids the essential singularity at x1:

Check:

Increasing the number of terms may improve accuracy. Here is a somewhat inaccurate approximation:

Check:

Increasing the number of terms produces a more accurate answer:

Increasing the number of terms further can produce nonsense due to numerical instability:

Combining an increase in the number of terms with a higher working precision often will reduce the error:

High-order derivatives with Method->EulerSum experience significant subtractive cancellation:

Using a higher working precision and additional terms produces an accurate answer:

For this problem, Method->NIntegrate with default options produces a correct answer:

Higher-order derivatives will again experience numerical instability:

Increasing WorkingPrecision will improve the accuracy:

An alternative is to increase the radius of the contour of integration: