function[output] = roots(a, b, c)
    D = b^2 - 4*a*c;
    if D < 0  
        disp("two complex")
        output = [(- b - sqrt(D)) / 2*a, (- b + sqrt(D)) / 2*a];
    elseif D == 0
        disp("one double")
        output = - b / 2*a;
    else
        disp("two complex")
        output = [(- b - sqrt(D)) / 2*a, (- b + sqrt(D)) / 2*a];
    end
end