I am interested in numerical methods to solve for the line which is tangent to a curve at exactly two points. I need a numerical method because the curve does not have a tidy explicit equation. I am more interested in robustness than speed, but (of course) speed is desirable.
This is part of larger algorithm development project: the curve in question has a parameter which, as it is varied, will cause the points of tangency to approach one another. At a unique value of x, y, and the parameter, the tangent lines become degenerate (merge); below this value of the parameter, there exists a unique bitangent to the curve, and at or above this value of the parameter there is no bitangent to the curve.
I have looked at related entries such as https://mathematica.stackexchange.com/questions/110668/how-to-find-a-tangent-line-with-2-points-of-tangency-for-a-curve and How can you find the equation of the line that is tangent at two distinct points to the curve? the former had a numerical method (or two?), but were written in the mathematic language and I wasn't sure what was going on or how to replicate the result in another language (such as julia or python), while the second is specific to polynomials and therefore not particularly helpful.
I know that I could solve a system of two equations involving my curve and two equations involving the first derivative of my curve, but the system appears ill-conditioned and I would very much like to find out if there is a better way first.
EDIT: The function is obviously not as simple as a quartic polynomial, but a quartic polynomial makes a decent model: for a quartic polynomial with two deep minima (such as y = x^4 -0.5x^3 + 9x^2 +2x -20), the tangency points are near the minima, but as the minima move toward one another, the intervening maximum becomes more and more shallow until the minima merge into one (such as y = x^4 + 3x^3 + 1.25x^2 + x + 5). At this point, however, there is still an existing bitangent (separated not by a maximum but by an inflection point). When the minima have moved sufficiently close to one another, however, the bitangent ceases to exist (I am pretty sure y = x^4 + 6x^3 + 14x^2 + 16x + 8 is an example). I want to find a numerical algorithm that will reliably find the bitangent if it exists and fail in a predictable way when it doesn't so that I can step through values of the parameter and find the threshold value at which the bitangent ceases to exist. This has applications, for example, in locating the plait point in ternary and higher phase diagrams.