function z = zerosquadpolstab(alpha,beta)
% MATLAB function computing the zeros of a quadratic polynomial
% $\xi\to \xi^2+\alpha\xi+\beta$ by means of the familiar discriminant
% formula $\xi_{1,2} = \frac{1}{2}(-\alpha+\sqrt{\alpha^2-4\beta})$ and
% \emph{Vieta's theorem}. This is a \emph{stable} implementation.
D = alpha^2-4*beta; % discriminant
if (D < 0), z = []; 
else
    wD = sqrt(D); 
    % Use discriminant formula only for zero far away from $0$ in order to
    % \emph{avoid cancellation} (German: Ausl"oschung). For the other zeros
    % use Vieta's formula. 
    if (alpha >= 0)
        t = 0.5*(-alpha-wD);
        z = [t;beta/t];
    else
        t = 0.5*(-alpha+wD);
        z = [beta/t;t];
    end
end

