  function z = AllRoots(alpha)
% z is a column vector consisting of all the roots of the
% function f(x) = sin(alpha*x)-x on the interval [-1 1]

% Observation 1. f(0) = 0.

% Observation 2. if f(z) = 0 then f(-z) = 0. So just look for the positive 
% roots in  [0 1].

% Observation 3. The roots fall in between the zeros of 
% f'(x) = alpha*cos(alpha*x)-1. Thus the critical values
% of f are those x that satisfy cos(alpha*x) = 1/alpha.
if alpha<1
    % f' is never zero so z = 0 is the only root..
    z = 0;
else
    z = 0;
    % The first critical value...
    xL = acos(1/alpha)/alpha;
    % The next critical value..
    xR = 2*pi/alpha - xL;
    % Since xR-xL = 2(pi-acos(1/alpha))/alpha >0  [xL,xR] is
    % a bracketing interval.
    fL = f(xL,alpha);
    fR = f(xR,alpha);
    while xL<=1 && fL*fR<0
        zk =fzero(@(x) f(x,alpha),[xL,xR]);
        % Is the root in [0,1]?
        if zk<=1
           % Don't forget that if zk is a root then so is -zk...
           z = [-zk z zk];
        end
        % The next critical value beyond xR is xL+2pi/alpha..
        L = xL; xL = xR; xR = L+2*pi/alpha; 
        fL = fR;
        fR = f(xR,alpha);
    end
end


function y =f(x,alpha)
y = sin(alpha*x)-x;

% -1 for not using that f is odd
% -1 for failure when alpha very close to 1 (alpha = 1.001)
% -1 for mistake in output vector (sometimes zero omitted, sometimes
% negative roots missing, sometimes duplicates)
% -1 for not handling alpha <= 1
