How to improve my curve fitting for an LED in spice
You mentioned professional models with a few components inside. A full device model certainly may include some AC components to deal with charge storage and inductance. But you are developing just the DC model. So don't worry about those details here.

The DC equation is a modified Shockley equation that looks something like this:

$$V=\eta\, V_T \operatorname{ln}\left(\cfrac{I}{I_S}+1\right)+I\cdot R_S$$

You want to figure out \$\eta\$, \$I_S\$, and \$R_S\$.



As \$V_T=\frac{k T}{q}\$ and is temperature-dependent, you do not want to complicate your life by making temperature yet another variable. So make sure that you control the diode's die temperature and keep it a constant. The best way to do this is to:


Allow the LED temperature time to equilibrate with the ambient temperature -- so if you are holding the LED for a while, give it time to cool down to room temperature, for example,
avoid making measurements over long periods of time where the ambient temperature has time to move,
make your measurements quickly to avoid heating the LED die,
avoid using higher currents that are likely to cause significant heating of the LED's die and, when necessary (it will be), use only very short pulse currents with long delays between (but not so long as to violate (1) or (2) above.)


(Or you can use an elevated, stable temperature heater to keep the LED at some fixed temperature.)



Let's look at the overall resistance illustrated by the above equation. This is found by computing \$\frac{\textrm{d} V}{\textrm{d} I}\$ of a slightly simplified version (ignoring the +1 term for now):

$$R \approx \frac{\eta\, V_T}{I} + R_S $$

That shows us that there are two parts to the resistance. And resistance is the slope of the \$I\$ vs \$V\$ curve. This breaks the curve into two parts, with a transition region between them. The two regions are where \$I \ll \frac{\eta\, V_T}{R_S}\$ and where \$I \gg \frac{\eta\, V_T}{R_S}\$.

There is a third region which can't be seen from this equation, but is present in the fuller equation presented earlier above. This is where \$I \approx I_S\$. You will need to stay away from this area because it is near here that the ignored +1 term matters more and complicates your life.

So you need to make at least two different measurements in at least two different regions. One where \$I_S\ll I \ll\frac{\eta\, V_T}{R_S}\$ (I'll call this Region II) and the other where \$I \gg\frac{\eta\, V_T}{R_S}\$ (Region III.) Region I is where the +1 term matters. You should stay away from Region I, entirely.

Region II may not heat up the LED die that much (good.) But Region III likely will, so you probably need pulse measurements for that region.

For Region II, you are attempting to figure out the value of \$\eta\$ and \$I_S\$. If you plot Region II measurements on a \$y=\operatorname{ln}\left(I\right)\$ vs \$x=V\$ curve, you will find that they fall on a fairly linear line. If not, you aren't fully within the region. Region I will curve away from the line at very low currents and Region III will also curve away from the line at very high currents. So you need to find and then make your Region II measurements in the linear area between Region I and Region III.

(There are transitions between these three Regions, too. Stay away from those transition areas, as well.)

Once you have these Region II measurements, \$\eta\$ is then inversely proportional to the slope of the plotted line. \$b=\operatorname{ln}\left(I_S\right)\$ is the y-axis intercept. If you make more than two such measurements at different values of \$I\$, then a standard linear least squares fitting routine can compute the slope (\$m\$) and intercept (\$b\$) for you. The slope you get from this will be \$m=\frac{1}{\eta\, V_T}\$, so just solve that: \$\eta=\frac{1}{m\cdot V_T}\$. The intercept will be in terms of \$\operatorname{ln}\left(I\right)\$, so you need to raise \$e\$ to that power to get \$I_S\$. Such fitting algorithms are everywhere. So you shouldn't have any trouble with that. (If you only make two measurements you can use simple algebra to figure out the resulting slope and intercept.)

To get \$R_S\$, you need to use large, pulsed currents in Region III, with those long delays between. Just two such measurements should be enough. (Again, to figure out where Region III is, you need to take a variety of pulsed measurements to probe that region and find out where it lays. Pick currents that require voltages well past the knee, here.)

Again, take at least two good measurements in Region III. More is okay, but then you are back to that least squares fitting routine again. From these you can figure out the slope. And the slope is the value of \$R_S\$.



That process should get you what you need. It uses relatively simple math tools (simple algebra for two measurements in each region, or linear least square fitting if more than two measurements in each region), too.

One thing I can see from your curve is that I don't think you have enough data to figure out \$R_S\$ well. You'll likely need higher current measurements.



To make the above concrete, imagine that \$V_T=26\:\textrm{mV}\$, \$\eta=3\$, \$I_S=1\:\textrm{pA}\$, and \$R_S=2\:\Omega\$. Region II is then \$1\:\textrm{pA}\ll I \ll 39\:\textrm{mA}\$ and Region III is \$I \gg 39\:\textrm{mA}\$.

Probe Region II at \$I_1=500\:\mu\textrm{A}\$ and \$I_2=3\:\textrm{mA}\$ (I picked those to create some interesting errors) to get \$V_1\approx 1.563\:\textrm{V}\$ and \$V_2\approx 1.708\:\textrm{V}\$. (I'm assuming we can only measure to four places.) On the log-lin plot, the slope is \$m=12.3569619\$ (It is valid to retain fuller precision for m [and for b, later].) Therefore, we compute \$\eta=\frac{1}{m V_T}\approx 3.1\$ -- which is close given the mV precision of the voltage measurement. The intercept is \$b=-26.9148338\$ and this leads to \$I_S=e^b\approx 2.05\:\textrm{pA}\$. This appears to be a more substantial error. But it is due to the effect of \$R_S\$, which at \$3\:\textrm{mA}\$ drops about \$6\:\textrm{mV}\$. Had I chosen a smaller high end current, such as \$1\:\textrm{mA}\$, I would have gotten about \$1.4\:\textrm{pA}\$ instead, which is closer.

You could then probe Region III with \$I_3=60\:\textrm{mA}\$ and \$I_4=100\:\textrm{mA}\$ (for very short periods) to get \$V_3=2.056\:\textrm{V}\$ and \$V_3=2.176\:\textrm{V}\$. (Again, these two currents are chosen close to one side of Region III to create interesting error, as well as to be practical instead of impractical.) From that, it is easy to find \$R_S\approx 3\:\Omega\$. Here, it's also wrong but not too far off. But you also need to take note that the dynamic part of the resistance was significant at \$1.34\:\Omega\$ and \$0.81\:\Omega\$, respectively. (Assuming our Region II parameter estimates.) So, of course there's an error. You could use that prior "knowledge" from Region II, namely our belief that \$\eta\approx 3.1\$ and \$I_S\approx 2.05\:\textrm{pA}\$, to then subtract about \$1\:\Omega\$ (a rough average of the two computed dynamic resistance values) to find \$R_S=2\:\Omega\$, which is of course much better.



All this should give you some ideas about how to figure things out and it manages to avoid overly complex mathematics by focusing on the dominant issues at hand in each region.