You can't do so. Well, not perfectly. Why not? Because I can always provide a function that has roots that are impossible to find. So I can sneak a pair of roots into any function that a rootfinder literally cannot find. It is not even that difficult to do so. The function will look quite normal and smooth, even a plot will not show anything strange, but if you knew exactly where to look, and to look very carefully, it would have an extra pair of roots. The trick is to start with a perfectly normal looking problem like the one you show, then add a term that is numerically close to a delta function, so narrow that you cannot see it on any plot.
That means your goal is technically impossible to solve using any numerical root finding method, even if a plot is used.
Next, I can give you a simple example where there are infinitely many roots in any small interval near zero. A simple example of that is the trivial function sin(1/x). You can NEVER write any expression for the first n roots above zero for that function, even though the roots are easy to write. If you pick any set of n roots for that function, I can trivially show there are infinitely many roots that are smaller than the ones you have provided.
Finally, there is the trivial issue - that you can look all you want for the first 3 roots of the function y(x)=x^2-1, but it only has two roots. And only one of them is positive.
The above reasons are why your problem is technically unsolvable. At the same time, it is practically solvable on most non-degenerate functions, to some extent. How?
Write code that evaluates your function at a sequence of points, starting from zero. On any interval, IF there is a zero crossing between consecutive pairs of points, then use fzero on that interval. Once you find a root, then continue looking forther out. Keep going until you hit n roots.
As long as the above sequence of points are spaced finely enough and you go far enough out, you will eventualy catch n roots (for any finite value of n, as long as there are at least n roots.) The problem is if you sample too coarsely, then an interval might have TWO zero crossings. And then you will not know to look at that interval. So in order to do the above solution, you need to know how finely to sample. Too coarse, and you miss roots. Too fine, and you waste a lot of CPU time, evaluating the function at points where there is no zero crossing.
The above scheme will work, in practice for a numerical solution, as long as your problem is not too degenerate. Could there be a symbolic code that will work? Sadly, it is impossible to solve some even simply written equations for a symbolic solution.