How to take the second derivative of a complicated function using limits [closed]

Question

I have a very complicated function of one variable consisting of exponential and error functions combined in various ways. Mathematica can't take the derivative of this function in the usual way (using prime or D) because it makes use of the function Im. Following the discussion here, I decided to use the limit definition of the derivative instead of prime or D:

dn[\[Delta]_] = Assuming[Element[\[Delta], Reals], Limit[(Realn[(\[Delta] + \[CurlyEpsilon])] - Realn[\[Delta]])/\[CurlyEpsilon], \[CurlyEpsilon] -> 0]];

This works for the first derivative, but returns a function which is infinity everywhere for the second derivative:

ddn[\[Delta]_] = 
  Assuming[Element[\[Delta], Reals], 
   Limit[(Realn[(\[Delta] + \[Xi])] - 2*Realn[\[Delta]] + 
       Realn[(\[Delta] - \[Xi])])/\[Xi]^2, \[Xi] -> 0]];

at any point. While trying to debug this problem, I discovered that multiplying the numerator and the denominator in the definition of the first derivative by 2 causes it too to give a function which is infinity at all points:

dn[\[Delta]_] = Assuming[Element[\[Delta], Reals], Limit[(2 Realn[(\[Delta] + \[CurlyEpsilon])] - 2 Realn[\[Delta]])/(2 \[CurlyEpsilon]), \[CurlyEpsilon] -> 0]];

and which I am utterly unable to explain. I can do both derivatives numerically using a finite number of points, and the results look reasonably well-behaved (i.e finite and smooth) within the region where the function is varying (it approaches zero at both infinities). Besides, the function made by summing terms which are all infinitely differentiable, so I don't see why the analytical method doesn't work...