% quadroot2.m roots of quadratic equation
% a* x^2 + b* x + c = 0

v= [a b c];
d= max(v); 
astore= a; bstore = b; cstore= c;
a = a/d; b= b/d; c= c/d;

if (a < eps) x1 = -c/b; sprintf('x1 = %24.15e', x1)
elseif (b^2 - 4*a*c <0) sprintf('complex roots'); 
else 
    if (b > 0) 
       x1= (-2*c) / (b + sqrt(b^2- 4*a*c)) ; 
       x2= (- b - sqrt(b^2 - 4*a*c))/(2*a);   
    else %  b < = 0
       x1= (- b + sqrt(b^2 - 4*a*c))/(2*a);
       x2= (-2*c) / (b - sqrt(b^2- 4*a*c)) ; 
    end 
    sprintf('x1 = %24.15e', x1)
    sprintf('x2 = %24.15e', x2)
end 

a = astore; b= bstore; c = cstore;