Guv:

I don't know who developed the solution that I used in my code. I got it from Perry's Chemical Engineers' Handbook.

The code I posted only worked for the case when R<0, which was when there were 3 real unequal roots.

I worked on it this morning and added the functionality to handle the other cases; R = 0 (three real roots, of which at least two are equal), R > 0 (one real root and two conjugate complex roots).

I tested (not extensively) the latest version against mathCAD's polyroots function and it worked for few examples I tried for each of the three states of R.

If anyone is interested, I could post the new code that now returns all the roots (complex and real).