You need to understand it is IMPOSSIBLE to know the number of solutions to a general nonlinear equation, or system of them.
In this case, in the interval of interest, there appear to be 3 real roots. But singe the problem is nonlinear, fzero cannot find them all. See my first comment. Accept it as truth, and know that since fzero can only interrogate the function at a finite number of points, it cannot know how to find all roots in that interval.
A simple solution is to evaluate the cuntion at a finite set of points yourself. So if we have:
x = linspace(0.00000000000001,0.99999999999999,1000);
So there are 3 spots identified where the function changes sign on that sampling.
idx = find(diff(sign(fx)))
and now we can find the three roots using fzero.
fzero(f,[x(idx(1)),x(idx(1) + 1)])
fzero(f,[x(idx(2)),x(idx(2) + 1)])
fzero(f,[x(idx(3)),x(idx(3) + 1)])