roots of functions with cusps (combined waves)

Question

If I have some non-negative functions with shapes similar to the one in the picture, (imagine, for simplicity, that each function is the sum of two $sin$ waves taken in absolute value.)

In order to find possible roots of the function, what is the most suitable algorithm that I can use and implement in a programming language? It seems that those cusps can "deceive" both fixed-point iteration and the Newton method. I was wondering if there are some "fixes" or variants to the iteration method to deal with those cusps and find where a cusp touches the x axis.

Answer 1

Taking care to choose a suitable step and tolerance, Newton-Raphson method works:

Zeros[f_, xmin_, xmax_, dx_, tol_] := Module[
  {g = D[f, x], sol = {}, xf = xmin + dx, xi = xmin, xm},
  While[xf < xmax,
        If[(f /. {x -> xi}) (f /. {x -> xf}) < tol,
           xm = (xi + xf) / 2.;
           While[Abs[f/g /. {x -> xm}] > tol,
                 xm = xm - f/g /. {x -> xm}];
           sol = Join[sol, {xm}]];
        xi = xf; xf = xi + dx];
  Union[sol, SameTest -> (Norm[#1 - #2] < tol &)]]

f[x_] := RealAbs[RealAbs[Sin[x^2 + x]] - Cos[x^2 - x]]
sol = Zeros[f[x], -5, 5, 10^-4, 10^-6]

pts = Transpose[{sol, ConstantArray[0, Length[sol]]}];
Plot[f[x], {x, -5, 5}, Epilog -> {Red, PointSize[0.012], Point[pts]}]

{-4.60497, -4.43113, -3.92699, -3.86297, -3.19534,
-2.93928, -2.34474, -1.98166, -0.785398, 0.785398,
0.886227, 1.53499, 2.93928, 3.19534, 3.86297,
3.92699, 4.0612, 4.2502, 4.77248, 4.9343}

Obviously nothing miraculous, you've to be careful, but not impossible either.