Borrowing almost verbatim from a recent response about finding [extrema](http://mathematica.stackexchange.com/questions/5575/how-to-find-all-the-local-minima-maxima-in-a-range/5591#5591), here is a method that is useful when your function is differentiable and hence can be "tracked" by NDSolve.

 f[x_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]]
 
 In[191]:= zeros = 
 Reap[soln = 
 y[x] /. First[
 NDSolve[{y'[x] == Evaluate[D[f[x], x]], y[10] == (f[10])}, 
 y[x], {x, 10, 0}, 
 Method -> {"EventLocator", "Event" -> y[x], 
 "EventAction" :> Sow[{x, y[x]}]}]]][[2, 1]]
 
 
 During evaluation of In[191]:= NDSolve::mxst: Maximum number of 10000 steps reached at the point x == 1.5232626281716416`*^-124. >>
 
 Out[191]= {{9.39114, 
 8.98587*10^-16}, {6.32397, -3.53884*10^-16}, {3.03297, \
 -8.45169*10^-13}, {0.886605, -4.02456*10^-15}}
 
 Plot[f[x], {x, 0, 10}, 
 Epilog -> {PointSize[Medium], Red, Point[zeros]}]

![enter image description here][1]

If it were a trickier function, one might use Method->{"Projection",...} to enforce the condition that y[x] is really the same as f[x]. This method may be useful in situations (if you can find them) where we have one function in one variable, and Reduce either cannot handle it or takes a long time to do so.

 [1]: https://i.sstatic.net/vFe5w.png