function [root1,root2]=gquad_solver(a,b,c)

% Function solves quadratic equation 
% ax^2+bx+c=0
% using modified formulae which takes into account 
% sign of b and prevents catastrophic cancellation errors

    d=discrim(a,b,c);
  
 % root of largest magnitude  
    q=-(b+sign(b)*d)/(2*a);
    root1=q;
    
 % other root
    root2=c/(a*root1);

end %gquad_solver()

% Besides solving the quadratic equation, this function 
% also provides elementary example of subfunction usage

function d = discrim(a,b,c);
    d = sqrt(b*b - 4*a*c);
end %discrim ()

