(1) When the series resistance is "small enough" such that \$v_D \approx V_1 \$, the current approximately follows the exponential diode current equation.
(2) When the series resistance is "large enough" such that \$v_{R1} \approx V_1 \$, the current approximately follows the linear Ohm's law equation.
So, you should only expect the current to approximately follow the diode exponential equation when \$v_{R1} << v_D\$.
This implies that you'll only see something resembling the exponential response for \$ V_1 < 0.8V \$ or so.
For this circuit, the diode current is given by:
\$i_D = I_S \exp(\frac{V_1 - i_DR_1}{nV_T})\$
or
\$nV_T\ln(\frac{i_D}{I_S}) + i_DR_1 = V_1\$
If you stare at this awhile, you see that the series current is approximately linear when:
\$i_D >> \dfrac{nV_T}{R_1}\ln(\frac{i_D}{I_S})\$
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.
Best Answer
A diode is a strongly non-linear electronic part. Hence, it is absolutely necessary to discriminate between static and dynamic (differential) resistances. Therefore, it is common practice to use capital letters (R) for static resistances and small letters (r) for diff. resistances. (By the way: Unfortunately, this general rule is violated by the contribution as contained in "learningaboutelectronics.com" as referenced in the first comment).
The first part of the equation under consideration is derived from the voltage-current relation of the pn junction, which is an exponential one and can be found in each textbook. It is derived from the differential quotient (slope) of the function Id=f(Vd) and contains with the "temperature voltage" Vt=26mV the temperature dependence of the unit (26mV for normal ambient temperatures). The second part (rb) is nothing else than an ohmic part (rather small) which contains the ohmic influence of the material between the external pins and the pn junction. This part has a linearizing effect for the Id=f(Vd) characteristics and comes into play for relatively large currents Id.
The static resistance of a diode is simply the ratio R=Vd/Id.