On the numerical roots of transcendent functions

Writing:

{a, b, c, d, e, f, g} = {-2, 2, -2, 2, 0, 6, 5} / 6;

h[t_] := (a + b t + Cos[e + f t])^2 + (c + d t + Sin[e + f t])^2 - g^2

sol1 = NSolve[h[t] == 0, Reals, WorkingPrecision -> 10][[All, 1, 2]]


instantly we get:

{-0.1995415116, 0.6420342510, 3.052979574, 4.390122149}

which are all real roots of this function as graphically verifiable:

Plot[h[t], {t, -25, 25}, AxesLabel -> {"t", "h[t]"}]


The best I can do to emulate NSolve[] is to apply the Bolzano's theorem coupled with the Newton-Raphson method:

sol2 = {};

T = 5;
dt = 10^-1;
tol = 10^-10;

ti = -T;
tf = -T + dt;

While[tf < T,
If[h[ti] h[tf] < tol,
t0 = (ti + tf) / 2.;
While[Abs[h[t0] / h'[t0]] > tol,
t0 = t0 - h[t0] / h'[t0]
];
sol2 = Join[sol2, {t0}]
];
ti = tf;
tf = tf + dt
]

sol1 == sol2


True

but it's evident that all depends heavily on the choice of parameters 0 < dt < T.

The question is therefore how to improve this approach or what other methods I could apply automatically to have the same solutions of NSolve[], as far as possible of course, I understand well that Mathematica algorithms will always be better than my "toys".

I was very fascinated by the methods suggested in the comments by @Michael E2. While the "graphical methods" are linked to the Plot[] functions implemented internally to Mathematica, the "spectral methods" are potentially exportable in less advanced spreadsheets.

In particular, aware of my limitations and without any pretense, implementing the Chebyshev-proxy rootfinder algorithm in the following way, the idea is to generate the coefficients of a polynomial that in a given interval has the same real roots as the assigned transcendent function, easily computable, for example, with Aberth-Ehrlich method, then adjust with the evergreen Newton-Raphson method:

(*Roots search range and tolerance*)
{tmin, tmax} = {-5, 5};
tol = 10^-10;

n = 10;
While[True,

(*Calculate Chebyshev's nodes*)
chebnodes = ConstantArray[0, n + 1];
For[i = 1, i <= n + 1, i++,
ti = Cos[Pi (i - 1) / n] (tmax - tmin) / 2 + (tmax + tmin) / 2.;
chebnodes[[i]] = h[ti]
];

(*Calculate Chebyshev's coefficients*)
chebcoeff = ConstantArray[0, n + 1];
For[i = 1, i <= n + 1, i++,
sum = 0;
For[j = 2, j <= n, j++,
sum = sum + Cos[Pi (i - 1) (j - 1) / n] chebnodes[[j]];
];
chebcoeff[[i]] = (chebnodes[[1]] + 2 sum + (-1)^(i-1) chebnodes[[n+1]])/n
];
chebcoeff[[1]] = chebcoeff[[1]] / 2;
chebcoeff[[n + 1]] = chebcoeff[[n + 1]] / 2;

(*Calculate error*)
err = 0;
For[i = 2 + n/2, i <= n + 1, i++,
err = err + Abs[chebcoeff[[i]]]
];

(*Check if it's necessary to double the nodes*)
If[err > tol, n = 2 n, Break[]]
];

(*Calculate coefficients of resolving polynomial*)
polycoeff = ConstantArray[0, n + 1];
For[i = 1, i <= n + 1, i++,
k = i;
For[j = 1, j <= (5 + 2 (n - i) - (-1)^(n + i)) / 4, j++,
If[i == 1,
l = 1,
l = 2^i (i + 2 j - 3) (i + j - 3)! / (4 (i - 1)! (j - 1)!)
];
polycoeff[[i]] = polycoeff[[i]] - (-1)^(n + j) l chebcoeff[[k]];
k = k + 2
]
];

(*Calculate polynomial's roots with Aberth-Ehrlich method*)
poly = 0;
For[i = 1, i <= n + 1, i++,
poly = poly + polycoeff[[i]] t^(i - 1)
];
polyroots = {ToRules[NRoots[poly == 0, t, Method -> "Aberth"]]}[[All, 1, 2]];

(*Select real roots within the interval*)
sol = ConstantArray[0, Length[polyroots]];
nsol = 0;
For[i = 1, i <= Length[polyroots], i++,
If[-1.001 < Re[polyroots[[i]]] < 1.001 && Abs[Im[polyroots[[i]]]] < 10^-15,
nsol = nsol + 1;
sol[[nsol]] = polyroots[[i]] (tmax - tmin) / 2 + (tmax + tmin) / 2
]
];

sol3 = ConstantArray[0, nsol];
For[i = 1, i <= nsol, i++,
t0 = sol[[i]];
While[Abs[h[t0] / h'[t0]] > tol,
t0 = t0 - h[t0] / h'[t0]
];
sol3[[i]] = t0
];

sol1 == sol3


True

It's undoubtedly a "toy code", from which, however, I had (and maybe others will have in the future) interesting ideas, most likely new notions compared to the classic "lessons for engineers".

• Related: spectral methods such as chebRoots[] and graphical methods such as FindAllCrossings[]. – Michael E2 Jun 3 at 15:08
• Is this question asking how NSolve works? Or do you have a specific use-case where it doesn't work, and you'd like help developing an approach for that? As far as I'm concerned, if NSolve can do it usually it means that anything you try won't be as good, so why even bother? – Fidel I. Schaposnik Jun 3 at 16:29
• @MichaelE2: I thank you infinitely, I have finally been able to acquire new information, although I'm aware that it will take time and a lot of study before I fully understand them! I have updated the request by attaching what I have learned again. – TeM Jun 4 at 18:06
• @FidelI.Schaposnik: Yes, I essentially asked for some robust strategy that NSolve[] could refer to. The pages linked by MichaelE2 were just what I was looking for. Of course, if you have any further information, they are welcome! – TeM Jun 4 at 18:07

NSolve for a given range of t solves your problem:
 NSolve[{h[t] == 0, -10 < t < 10}, t]