Too late, since you have an answer that solves your problem to your satisfaction. I want to point out something you probably need to know though.
Finding the first n solutions using a numerical solver is actually an impossible task to solve. This is because you cannot insure that you have found all solutions over some interval. You may lose solutions, and that means you got close to the first n solutions, but you missed a few in there. And sometimes, a numerical solver can find solutions that are not in fact solutions. The problem is fraught with issues that touch on the weak points of numercal solvers, and if you throw a difficult problem at them, well, OOPS!
Ugh. Why am I saying this? First, how does a tool like fzero work? I'll suggest fzero is your best choice for a 1-d numerical rootfinder, since it is more robust than a tool like fsolve. Fzero prefers a bracket, where the function is known to cross zero. So it wants an interval [a,b], such that f(a) and f(b) have different signs. (Or f(x) could be exactly zero at one of the endpoints. That is ok too.) If you don't give fzero a bracket, but only ONE point (a), then fzero will search for a second point where f(x) has a different sign than f(a). And then once it has a bracket, it is happy again.
But once fzero has a bracket, then it uses a hybrid set of methods that insure efficient convergence on good problems, but also insure the convergence cannot go too slowly on some classically known bad problems that can arise. We don't need to worry about those tweaks, but only accept that a tool like fsolve does not use them, and so can have issues of its own. (A good course on numerical analysis would be nice here, but I won't be spending the time to teach it. Sorry.)
Let me suggest a simple problem to solve, then I'll modify it in just a tiny way. Find the root of the function
Yes, I know the root lies at x==100. But the computer does not. Numerical tools do NOT actually look at the function itself, but only as if the function is a black box they cannot look inside. So they send in x values, and then look at what comes out of the magical black box.
opts = optimset('fzero');
[x,fval,exitflag] = fzero(f1,x0,opts)
Search for an interval around 1 containing a sign change:
Func-count a f(a) b f(b) Procedure
1 1 0.99 1 0.99 initial interval
3 0.971716 0.990283 1.02828 0.989717 search
5 0.96 0.9904 1.04 0.9896 search
7 0.943431 0.990566 1.05657 0.989434 search
9 0.92 0.9908 1.08 0.9892 search
11 0.886863 0.991131 1.11314 0.988869 search
13 0.84 0.9916 1.16 0.9884 search
15 0.773726 0.992263 1.22627 0.987737 search
17 0.68 0.9932 1.32 0.9868 search
19 0.547452 0.994525 1.45255 0.985475 search
21 0.36 0.9964 1.64 0.9836 search
23 0.0949033 0.999051 1.9051 0.980949 search
25 -0.28 1.0028 2.28 0.9772 search
27 -0.810193 1.0081 2.81019 0.971898 search
29 -1.56 1.0156 3.56 0.9644 search
31 -2.62039 1.0262 4.62039 0.953796 search
33 -4.12 1.0412 6.12 0.9388 search
35 -6.24077 1.06241 8.24077 0.917592 search
37 -9.24 1.0924 11.24 0.8876 search
39 -13.4815 1.13482 15.4815 0.845185 search
41 -19.48 1.1948 21.48 0.7852 search
43 -27.9631 1.27963 29.9631 0.700369 search
45 -39.96 1.3996 41.96 0.5804 search
47 -56.9262 1.56926 58.9262 0.410738 search
49 -80.92 1.8092 82.92 0.1708 search
51 -114.852 2.14852 116.852 -0.168524 search
Search for a zero in the interval [-114.852, 116.852]:
Func-count x f(x) Procedure
51 116.852 -0.168524 initial
52 100 0 interpolation
Zero found in the interval [-114.852, 116.852]
So fzero does eventialluy look out as far as x==100. once it gets there, it zooms right in on the solution.
But if we give fzero a bracket, then it is supremely happy, even if the bracket is very wide.
[x,fval,exitflag] = fzero(f1,x0,opts)
Func-count x f(x) Procedure
2 0 1 initial
3 100 1.11022e-16 interpolation
4 100 1.11022e-16 interpolation
Zero found in the interval [0, 1000]
But now, suppose we give fzero a slightly different function?
f2 = @(x) 1 - 0.01*x - exp(-100*(x - 23).^2);
f1 and f2 are virtually identical everywhere, except for a tiny blip, that occurs only around x==23.
So f2 actually has a zero crossing, TWO of tem, in fact, very close to x==23. And everywhere ealse, f1 and f2 are pretty much identical. I could have made the blip an even sharper one too. Visualize an approximate Dirac delta function.
The problem is, no numerical tool can find that pair of roots, especially if I make the blip a very sharp one. If I give fzero the function f2, it will still find only the root at x==100.
[x,fval,exitflag] = fzero(f2,x0,opts)
Search for an interval around 1 containing a sign change:
Func-count a f(a) b f(b) Procedure
1 1 0.99 1 0.99 initial interval
3 0.971716 0.990283 1.02828 0.989717 search
5 0.96 0.9904 1.04 0.9896 search
7 0.943431 0.990566 1.05657 0.989434 search
9 0.92 0.9908 1.08 0.9892 search
11 0.886863 0.991131 1.11314 0.988869 search
13 0.84 0.9916 1.16 0.9884 search
15 0.773726 0.992263 1.22627 0.987737 search
17 0.68 0.9932 1.32 0.9868 search
19 0.547452 0.994525 1.45255 0.985475 search
21 0.36 0.9964 1.64 0.9836 search
23 0.0949033 0.999051 1.9051 0.980949 search
25 -0.28 1.0028 2.28 0.9772 search
27 -0.810193 1.0081 2.81019 0.971898 search
29 -1.56 1.0156 3.56 0.9644 search
31 -2.62039 1.0262 4.62039 0.953796 search
33 -4.12 1.0412 6.12 0.9388 search
35 -6.24077 1.06241 8.24077 0.917592 search
37 -9.24 1.0924 11.24 0.8876 search
39 -13.4815 1.13482 15.4815 0.845185 search
41 -19.48 1.1948 21.48 0.7852 search
43 -27.9631 1.27963 29.9631 0.700369 search
45 -39.96 1.3996 41.96 0.5804 search
47 -56.9262 1.56926 58.9262 0.410738 search
49 -80.92 1.8092 82.92 0.1708 search
51 -114.852 2.14852 116.852 -0.168524 search
Search for a zero in the interval [-114.852, 116.852]:
Func-count x f(x) Procedure
51 116.852 -0.168524 initial
52 100 0 interpolation
Zero found in the interval [-114.852, 116.852]
In fact, it evaluated the function at x=21.48 and at x=29.9631. It completely missed the pair of roots near x=23.
And, fzero can also be tricked, if you give it a function with a zero crossing that is not actually a root.
[x,fval,exitflag] = fzero(f3,x0,opts)
Func-count x f(x) Procedure
2 1 1.55741 initial
3 2 -2.18504 interpolation
4 1.85165 -3.46644 interpolation
5 1.41615 6.4146 bisection
6 1.47895 10.8572 interpolation
7 1.6339 -15.8261 bisection
8 1.61954 -20.4982 interpolation
9 1.58798 -58.1765 bisection
10 1.55642 69.5775 bisection
11 1.55783 77.1347 interpolation
12 1.56502 173.074 bisection
13 1.56861 457.679 bisection
14 1.5722 -710.251 bisection
15 1.57182 -980.91 interpolation
16 1.57041 2574.06 bisection
17 1.57111 -3169.73 interpolation
18 1.57104 -4124.11 interpolation
19 1.57088 -11801.6 bisection
20 1.57072 13697.2 bisection
21 1.57073 14893.3 interpolation
22 1.57077 32636.8 bisection
23 1.57078 80720.7 bisection
24 1.5708 -170547 bisection
25 1.57079 306518 interpolation
26 1.5708 -384460 interpolation
27 1.5708 -515558 interpolation
28 1.5708 1.51194e+06 bisection
29 1.5708 -1.56465e+06 interpolation
30 1.5708 -1.62117e+06 interpolation
31 1.5708 -3.36385e+06 bisection
32 1.5708 -7.27282e+06 bisection
33 1.5708 -1.73587e+07 bisection
34 1.5708 4.48791e+07 bisection
35 1.5708 -5.66155e+07 interpolation
36 1.5708 -7.66641e+07 interpolation
37 1.5708 2.16494e+08 bisection
38 1.5708 -2.37393e+08 interpolation
39 1.5708 -2.62758e+08 interpolation
40 1.5708 -5.88385e+08 bisection
41 1.5708 -1.54689e+09 bisection
42 1.5708 2.45914e+09 bisection
43 1.5708 3.48749e+09 interpolation
44 1.5708 -8.33979e+09 bisection
45 1.5708 1.19881e+10 interpolation
46 1.5708 2.13107e+10 interpolation
47 1.5708 -2.74039e+10 bisection
48 1.5708 -3.1975e+10 interpolation
49 1.5708 -7.67527e+10 bisection
50 1.5708 1.91689e+11 bisection
51 1.5708 -2.56007e+11 interpolation
52 1.5708 -3.85261e+11 interpolation
53 1.5708 7.63028e+11 bisection
54 1.5708 1.49707e+12 interpolation
55 1.5708 -1.55633e+12 bisection
56 1.5708 -1.58761e+12 interpolation
57 1.5708 -3.24064e+12 bisection
58 1.5708 -6.76496e+12 bisection
59 1.5708 -1.48279e+13 bisection
60 1.5708 -3.64003e+13 bisection
61 1.5708 7.86301e+13 bisection
62 1.5708 -1.33542e+14 interpolation
63 1.5708 1.93489e+14 interpolation
64 1.5708 3.66869e+14 interpolation
65 1.5708 -4.19946e+14 bisection
66 1.5708 -5.83048e+14 interpolation
67 1.5708 -1.20927e+15 bisection
Current point x may be near a singular point. The interval [1, 2]
reduced to the requested tolerance and the function changes sign in the interval,
but f(x) increased in magnitude as the interval reduced.
So fzero finds a root at pi/2, but then it returns an exitfkag of -5, to say that it thinks this is not actlually a root. An exitflag of -5 means it thinks this is a point of singularity, but not a root.
help fzero
FZERO Single-variable nonlinear zero finding.
X = FZERO(FUN,X0) tries to find a zero of the function FUN near X0,
if X0 is a scalar. It first finds an interval containing X0 where the
function values of the interval endpoints differ in sign, then searches
that interval for a zero. FUN is a function handle. FUN accepts real
scalar input X and returns a real scalar function value F, evaluated
at X. The value X returned by FZERO is near a point where FUN changes
sign (if FUN is continuous), or NaN if the search fails.
X = FZERO(FUN,X0), where X0 is a vector of length 2, assumes X0 is a
finite interval where the sign of FUN(X0(1)) differs from the sign of
FUN(X0(2)). An error occurs if this is not true. Calling FZERO with a
finite interval guarantees FZERO will return a value near a point where
FUN changes sign.
X = FZERO(FUN,X0), where X0 is a scalar value, uses X0 as a starting
guess. FZERO looks for an interval containing a sign change for FUN and
containing X0. If no such interval is found, NaN is returned.
In this case, the search terminates when the search interval
is expanded until an Inf, NaN, or complex value is found. Note: if
the option FunValCheck is 'on', then an error will occur if an NaN or
complex value is found.
X = FZERO(FUN,X0,OPTIONS) solves the equation with the default optimization
parameters replaced by values in the structure OPTIONS, an argument
created with the OPTIMSET function. See OPTIMSET for details. Used
options are Display, TolX, FunValCheck, OutputFcn, and PlotFcns.
X = FZERO(PROBLEM) finds the zero of a function defined in PROBLEM.
PROBLEM is a structure with the function FUN in PROBLEM.objective,
the start point in PROBLEM.x0, the options structure in PROBLEM.options,
and solver name 'fzero' in PROBLEM.solver.
[X,FVAL]= FZERO(FUN,...) returns the value of the function described
in FUN, at X.
[X,FVAL,EXITFLAG] = FZERO(...) returns an EXITFLAG that describes the
exit condition. Possible values of EXITFLAG and the corresponding exit
conditions are
1 FZERO found a zero X.
-1 Algorithm terminated by output function.
-3 NaN or Inf function value encountered during search for an interval
containing a sign change.
-4 Complex function value encountered during search for an interval
containing a sign change.
-5 FZERO may have converged to a singular point.
-6 FZERO can not detect a change in sign of the function.
[X,FVAL,EXITFLAG,OUTPUT] = FZERO(...) returns a structure OUTPUT
with the number of function evaluations in OUTPUT.funcCount, the
algorithm name in OUTPUT.algorithm, the number of iterations to
find an interval (if needed) in OUTPUT.intervaliterations, the
number of zero-finding iterations in OUTPUT.iterations, and the
exit message in OUTPUT.message.
Examples
FUN can be specified using @:
X = fzero(@sin,3)
returns pi.
X = fzero(@sin,3,optimset('Display','iter'))
returns pi, uses the default tolerance and displays iteration information.
FUN can be an anonymous function:
X = fzero(@(x) sin(3*x),2)
FUN can be a parameterized function. Use an anonymous function to
capture the problem-dependent parameters:
myfun = @(x,c) cos(c*x); % The parameterized function.
c = 2; % The parameter.
X = fzero(@(x) myfun(x,c),0.1)
Limitations
X = fzero(@(x) abs(x)+1, 1)
returns NaN since this function does not change sign anywhere on the
real axis (and does not have a zero as well).
X = fzero(@tan,2)
returns X near 1.5708 because the discontinuity of this function near the
point X gives the appearance (numerically) that the function changes sign at X.
See also ROOTS, FMINBND, FUNCTION_HANDLE.
Documentation for fzero
doc fzero
Other uses of fzero
optim/fzero
These are general problems of such a tool. We can always trick a solver, as long as we understand how the tool works. So insuring you have found the first n roots, this is actually an impossible task, as I said before.
It gets worse of course. Consider the function f4.
For the function f3(x), there are infinitely many roots in the interval [0,b], for ANY positive value of b, no matter how small you make b. So again, finding the first n positive roots of f3 is a mathematical impossibility.
There are things you can do to try to make the problem semi-tractable. I posted a function called allcrossings on the file exchange some time ago, but allcrossings is not perfect. allcrossings uses a scheme where it evaluates the function at many points on an interval, looking for any sign change. Then on any sub-interval, it calls fzero, because it has a bracket that should contain a root. But allcrossings does not find the first n roots. It tries to find ALL roots, within limits on its abilities.
The point is? Finding the first n roots is a difficult problem for a numerical solver, not one for which there will ever be an easy solution.
