Diode equation with a series resistance
I'm going to assume you are asking about a diode without any external circuit other than an applied voltage, and that you are asking specifically about how to handle the non-zero simulation parameter, \$R_S\$, in computing the total diode current given an impressed diode voltage.

In that case, you have the following puzzle:

$$I_D=I_S\cdot\left(e^{\frac{V_D-R_S\cdot I_D}{n\cdot V_T}}-1\right)$$

You can solve this, iteratively. For example, suppose the applied voltage is \$V_D=600\:\textrm{mV}\$, the thermal voltage is \$V_T=26\:\textrm{mV}\$, the parasitic resistance is \$R_S=600\:\textrm{m}\Omega\$, the emission coefficient is \$n=1.7\$, and the saturation current is \$I_S=2.5\:\textrm{nA}\$. Then, you can guess at the current by first ignoring (setting it to zero) the term in the exponential, namely \$R_S\cdot I_D\$, and then iteratively compute it -- substituting in the new value each time. If you do this, you will get \$I_D=1.9144894\:\textrm{mA}\$ as a result.

However, you can instead decide to drop the \$-1\$ term in the standard equation and solve it as follows:

$$V_D-R_S\cdot I_D \approx n\cdot V_T\cdot\operatorname{ln}\left(\frac{I_D}{I_S}\right)$$

The above is approximate, because we've neglected the \$-1\$ term. But it is negligible in most cases since the exponential (for all practical currents) dominates the results. Besides, these simulation parameters aren't physical constants of the universe. They are just simulation parameters. So the above equation is almost always sufficiently close.

If you solve that equation for \$I_D\$ you can use the Lambert-W function and this results in:

$$I_D = \frac{n\cdot V_T}{R_s}\cdot\operatorname{LambertW}\left[\frac{I_S\cdot R_S}{n\cdot V_T}\cdot e^{\frac{V_D}{n\cdot V_T}}\right]$$

(If you are interested in what the LambertW function is [how it is defined] and in seeing a fully worked example on how to apply it to solving problems like these, then see: Differential and Multistage Amplifiers(BJT).)

Plugging in the values, I get \$I_D=1.914491841\:\textrm{mA}\$. Which is, of course, very very close to the earlier iterative solution.

Note also that dimensional analysis works out in the above equations. As it should.



Bottom line is that you can either iterate or else you can develop a closed solution using the Lambert-W function (also known as the product-log function.) Both methods work reasonably well.