Somewhat similar to my previous wonderings on TI's CORDIC implementation, today I'm wondering if it's known precisely what algorithm the numeric solver (MATH -> Solver...) on the 83+ family uses.

I have a vague memory of somebody claiming it used Newton's method at some point, but I think it's more likely to use the bisection method, since EOS is unable to symbolically construct the derivative of a function to find a root of. The secant method is a feasible alternative, but seems suboptimal because it successfully converging on a solution depends strongly on the initial estimate (and it requires two points to estimate at).

That the solver can sometimes generate ERR: NO SIGN CHANGE seems like a point in favor of it using the bisection method, since "no sign change" seems to imply it needs to find two points on the specified function that have opposite sign which is also a requirement of the bisection method. If that's true however, it's unclear how it goes about finding those points initially. For many functions it could look at the signs of the initial guess and upper/lower bounds to search for a zero within, but it still seems to have no trouble with functions that have the same sign at all three of those points: for instance -(x+6)²+4 is solved successfully with the default range (-1e99 to 1e99) and a guess of 10, where the function has a negative value at all three of those X values.

Does anybody have ideas or primary sources regarding how the solver is implemented? I'd be interested to hear!

Continuing my literature survey beyond Wikipedia, the ITP method seems interesting, but it's much newer than the 8x calculators and still seems dependent on finding a sign change. The question remains of how EOS locates a sign change, which it seems to need.

Turning to Numerical Methods for Engineers, chapters 5 and 6 offer nicely-presented suggestions. Chapter 5 includes the bisection and false-position methods, which as discussed first require finding a range of a target function that has opposite sign at either end. Going into "open methods", we find the algorithms that don't require an initial known range containing a sign change. This includes the secant method and Newton's method as we've already looked at, but also suggests Brent's method that seems promising. MATLAB's documentation also notes that its implementation won't find zeroes that don't cross the X axis, and I'm not sure offhand if TI's will- if it does, that could indicate it does something different.

Brent's method is documented as the method MATLAB uses in its fzero built-in function, which seems promising and it also turns out that TI|BD also claims that solve( uses Brent's method. It's not clear where the initial guess require by EOS would be used here however, unless it's used as the initial midpoint in inverse quadratic interpolation: in general, Brent's method doesn't use an initial guess so there's some difference between the typical formulation and what TI does if indeed EOS uses Brent's method.

I think it's possible that TI uses "no sign change" as a teacher-and-student understandable name for a different failure mode of the solver, even if their solver doesn't actually require a sign change.

The TI-Nspire and other calculators use a "combination of bisection and the secant method", whatever that means.
I consider it likely TI reused this algorithm. I've since purchased the book for all of three dollars to fuel my interest in 1970s computing and engineering; I'll let you know what I find.

Another algorithm is described for the 85/86, but I haven't evaluated if it is the same algorithm or a different one. Either way, TI seems to be willing to give out this information- consider asking?

The 85/86 approach you link sounds like it may be the same as the one mentioned for the Nspire, with the 85/86 one being described in more detail.

The discussion of finding a sign change is particularly interesting:This doesn't very clearly specify the heuristics, but does suggest that the "no sign change" error indicates the heuristic search for a range containing a zero failed.

iPhoenix has been able to share images of the copy of Shampine and Allen's book he got, which is moderately interesting. I haven't spent much time studying the implementation, but it seems plausible that their ZEROIN function is exactly the algorithm used by TI.

In the prose describing this implementation, they say that in general the two endpoints B and C must bracket a zero (f(B)*f(C)<0) which is as we'd expect from the normal closed methods, consistent with the description of this algorithm as combining bisection and the secant method. They go on to describe "abnormal" situations though, which include cases where the endpoints do not have opposite sign:(Recall that a root of even multiplicity is a point x where f(x)=0, but all points on either side of x have the same sign: the function does not cross the X axis.)
Here's the code I manually transcribed (because OCR on this FORTRAN didn't work very well). It's possible there are still some errors, because I don't fully understand how to read archaic FORTRAN like this.
Code:
This seems like it ticks the major boxes: it still works if the provided endpoints have the same sign, and gives up after a while if unable to find a sign change (ERR: NO SIGN CHANGE). It doesn't permit the caller to provide an initial guess, but the TI-83 Plus manual appears to provide us with enough information to understand how it gets incorporated:
So it seems like the initial guess would replace the midpoint between the two endpoints as the point at which to bisect; the first iteration won't strictly be bisecting, depending on what guess is provided.

A moderately interesting prediction made by these observations is that because roots of even multiplicity cannot be exactly pinpointed, it should be relatively easy to construct functions where the calculator returns answers that are off by a few least-significant digits, depending on the permitted relative and absolute errors.

Just verifying that I kind of understand how to compile and run the FORTRAN code we got, add the following main function (the example for the function's use provided in the book) to a file zeroin.f then compile as gfortran zeroin.f.

Code:
Running it on my x86 machine prints a slightly different result from what the book says, but it's an accurate result. I think the difference is that the book expects things to be run on an IBM 360/67 that appears to have used 32-bit HFP format, whereas I'll be using an IEEE 754 single-precision float today (or possibly double-precision, but I think the type for an implicit float is still REAL(4) which should be single-precision). My machine prints:
Code:
whereas the book says the result printed is
Code:
where the difference is that the IBM machine finds what it considers to be an exact result, even though both are off the infinite-precision exact result by around 3e-7.