Solving diode circuit with iteration - why doesn't it work?

Your electrical analysis is fine, it is really just a question of which numerical method to use, since the equation cannot be solved analytically.

Different numerical methods are appropriate for different situations. I'm not quite sure what the problem with using your method is in this circuit, but I think it is because the slope is so steep with the 10k resistor that the next estimate overshoots the correct value by way too much.

A different numerical method that can work is basically a binary search. The key to notice is that the function \$A - B e^{x}\$ is strictly decreasing, so that if the left hand \$V_D\$ is too low, you know the right hand \$V_D\$ is too high. By hand the process looks something like this (Python):

>>> def vd1(vd0):
...   return 5-6.9144e-12*math.exp(vd0/0.025)
...
>>> vd1(0.7)
-4.999999845336939
>>> vd1(0.6)
4.8168436139454105
>>> vd1(0.65)
3.646647188565237
>>> vd1(0.675)
1.3212056451829235
>>> vd1(0.68)
0.5067104283223642
>>> vd1(0.678)
0.8521709473357619
>>> vd1(0.679)
0.6828948324788575
>>> vd1(0.6791)
0.6655918288722198
>>> vd1(0.67905)
0.6742519821744661
>>> vd1(0.67902)
0.6794397665027283

So \$V_D \approx 0.679\$V.