Very recently, I learned a useful procedure due to J. P. Boyd (see also this and this) that involves expanding a function as a Chebyshev polynomial series, forming the so-called "colleague matrix", and then finding the eigenvalues of this matrix, which are often good approximations to the roots of the original function. I shall present how to do a barebones version of this strategy in Mathematica.

The first step is to form the Chebyshev series coefficients of the function concerned. One could either start from the Taylor series and use any number of methods to convert these coefficients to Chebyshev coefficients. Another route, which I shall present here, is to evaluate the function at the so-called "Chebyshev nodes", suitably transformed to the domain of the original function, and then using the (type-1) discrete cosine transform to produce the Chebyshev coefficients. It goes like this:

f[x_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]]; {xmin, xmax} = {1/2, 100}; n = 128; (* arbitrarily chosen integer *) prec = 20; (* precision *) cnodes = Rescale[N[Cos[Pi Range[0, n]/n], prec], {-1, 1}, {xmin, xmax}]; cc = Sqrt[2/n] FourierDCT[f[cnodes], 1]; cc[[{1, -1}]] /= 2; 

n here should be chosen to be greater than or equal to the number of roots known to be in the given interval. That is usually determined through a preliminary Plot[]. (A more sophisticated version, as used in the Chebfun package for MATLAB, does an adaptive increase in the order of the underlying Chebyshev approximation using a technique similar to Clenshaw-Curtis quadrature. I have chosen not to implement this for now.)

With the Chebyshev coefficients now available, one can then form the colleague matrix like so:

colleague = SparseArray[{{i_, j_} /; i + 1 == j :> 1/2, {i_, j_} /; i == j + 1 :> 1/(2 - Boole[j == 1])}, {n, n}] - SparseArray[{{i_, n} :> cc[[i]]/(2 cc[[n + 1]])}, {n, n}]; 

We can then find the eigenvalues of this matrix, filter out irrelevant eigenvalues, and then use Rescale[] to transform the remaining eigenvalues back to the domain of the original function:

rts = Sort[Select[DeleteCases[ Rescale[Eigenvalues[comrade], {-1, 1}, {xmin, xmax}], _Complex | _DirectedInfinity], xmin <= # <= xmax &]] {0.8866466, 3.0327390, 6.3243343, 9.3928065, 12.5892495, 15.6908929, 18.8634014, 21.9762069, 25.1562125, 28.3003285, 31.4743983, 34.5197235, 37.6746911, 40.8716527, 43.9572481, 47.1394116, 50.2776398, 53.3723848, 56.5916636, 59.6440741, 62.8638369, 66.0001767, 69.0778338, 72.2121395, 75.3874360, 78.5335813, 81.6858635, 84.8171664, 87.9688503, 91.1015749, 94.2504610, 97.3866141} 

The eigenvalues obtained look a bit rough here; this is due to the ill-behavior of the original function near the origin. For more well-behaved functions, the eigenvalues are often quite accurate roots of the original function. Of course, one can always use FindRoot[] to polish off these approximations.